home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1993 July / InfoMagic USENET CD-ROM July 1993.ISO / sources / misc / volume7 / ephem / part03 < prev    next >
Encoding:
Text File  |  1989-06-03  |  59.3 KB  |  2,146 lines
  1. Newsgroups: comp.sources.misc
  2. From: allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc)
  3. Subject: v07i011: astronomical ephemeris - 3 of 3
  4. Keywords: ephemeris astro
  5. Organization: Dimensional Medicine, Inc.  Minnetonka, MN.
  6. Reply-To: ecd@ncs-med.UUCP (Elwood C. Downey)
  7.  
  8. Posting-number: Volume 7, Issue 11
  9. Submitted-by: ecd@ncs-med.UUCP (Elwood C. Downey)
  10. Archive-name: ephem/part03
  11.  
  12. [BTW, not putting "guard characters" at the beginnings of lines in a shar'ed
  13. file is asking for trouble on some systems.  Especially with idiot UUCP mail
  14. systems -- some, perhaps all, of this group *does* get sent over mail links.
  15. One "From" will produce a mangled transmission.  And Bitnet seems to know of
  16. other ways to destroy such files.  ++bsa]
  17.  
  18. #!/bin/sh
  19. echo mkdir lib
  20. mkdir lib
  21. echo cd lib
  22. cd lib
  23. echo extracting Makefile
  24. cat > Makefile << 'xXx'
  25. CLNFLAGS=$(CLNF)
  26. LNFLAGS=$(CLNFLAGS) $(LNF)
  27. CFLAGS=$(CLNFLAGS) $(CF)
  28. LINTFLAGS=$(CFLAGS) $(LINTF)
  29. LINTLIBS=
  30.  
  31. .c.a:
  32.     -$(CC) -c $(CFLAGS) $<
  33.  
  34. .PRECIOUS:    lib.a
  35.  
  36. lib.a:    lib.a(aa_hadec.o) \
  37.     lib.a(anomaly.o) \
  38.     lib.a(cal_mjd.o) \
  39.     lib.a(eq_ecl.o) \
  40.     lib.a(moon.o) \
  41.     lib.a(moonnf.o) \
  42.     lib.a(moonrs.o) \
  43.     lib.a(nutation.o) \
  44.     lib.a(obliq.o) \
  45.     lib.a(parallax.o) \
  46.     lib.a(pelement.o) \
  47.     lib.a(plans.o) \
  48.     lib.a(precess.o) \
  49.     lib.a(refract.o) \
  50.     lib.a(riset.o) \
  51.     lib.a(sex_dec.o) \
  52.     lib.a(sun.o) \
  53.     lib.a(sunrs.o) \
  54.     lib.a(utc_gst.o)
  55.     ar rv $@ $?
  56.     ranlib $@
  57.     rm -f $?
  58. xXx
  59. echo extracting aa_hadec.c
  60. cat > aa_hadec.c << 'xXx'
  61. #include <stdio.h>
  62. #include <math.h>
  63. #include "astro.h"
  64.  
  65. /* given latitude (n+, radians), lat, altitude (up+, radians), alt, and
  66.  * azimuth (angle round to the east from north+, radians),
  67.  * return hour angle (radians), ha, and declination (radians), dec.
  68.  */
  69. aa_hadec (lat, alt, az, ha, dec)
  70. float lat;
  71. float alt, az;
  72. float *ha, *dec;
  73. {
  74.     aaha_aux (lat, az, alt, ha, dec);
  75. }
  76.  
  77. /* given latitude (n+, radians), lat, hour angle (radians), ha, and declination
  78.  * (radians), dec,
  79.  * return altitude (up+, radians), alt, and
  80.  * azimuth (angle round to the east from north+, radians),
  81.  */
  82. hadec_aa (lat, ha, dec, alt, az)
  83. float lat;
  84. float ha, dec;
  85. float *alt, *az;
  86. {
  87.     aaha_aux (lat, ha, dec, az, alt);
  88. }
  89.  
  90. /* the actual formula is the same for both transformation directions so
  91.  * do it here once for each way.
  92.  * N.B. all arguments are in radians.
  93.  */
  94. static
  95. aaha_aux (lat, x, y, p, q)
  96. float lat;
  97. float x, y;
  98. float *p, *q;
  99. {
  100.     static float lastlat = -1000.;
  101.     static float sinlastlat, coslastlat;
  102.     float sy, cy;
  103.     float sx, cx;
  104.     float sq, cq;
  105.     float a;
  106.     float cp;
  107.  
  108.     /* latitude doesn't change much, so try to reuse the sin and cos evals.
  109.      */
  110.     if (lat != lastlat) {
  111.         sinlastlat = sin (lat);
  112.         coslastlat = cos (lat);
  113.         lastlat = lat;
  114.     }
  115.  
  116.     sy = sin (y);
  117.     cy = cos (y);
  118.     sx = sin (x);
  119.     cx = cos (x);
  120.  
  121. /* define GOODATAN2 if atan2 returns full range -PI through +PI.
  122.  */
  123. #ifdef GOODATAN2
  124.     *q = asin ((sy*sinlastlat) + (cy*coslastlat*cx));
  125.     *p = atan2 (-cy*sx, -cy*cx*sinlastlat + sy*coslastlat);
  126. #else
  127. #define    EPS    (1e-20)
  128.     sq = (sy*sinlastlat) + (cy*coslastlat*cx);
  129.     *q = asin (sq);
  130.     cq = cos (*q);
  131.     a = coslastlat*cq;
  132.     if (a > -EPS && a < EPS)
  133.         a = a < 0 ? -EPS : EPS; /* avoid / 0 */
  134.     cp = (sy - (sinlastlat*sq))/a;
  135.     if (cp >= 1.0)    /* the /a can be slightly > 1 */
  136.         *p = 0.0;
  137.     else if (cp <= -1.0)
  138.         *p = PI;
  139.     else
  140.         *p = acos ((sy - (sinlastlat*sq))/a);
  141.     if (sx>0) *p = 2.0*PI - *p;
  142. #endif
  143. }
  144. xXx
  145. echo extracting anomaly.c
  146. cat > anomaly.c << 'xXx'
  147. #include <stdio.h>
  148. #include <math.h>
  149. #include "astro.h"
  150.  
  151. #define    TWOPI    (2*PI)
  152.  
  153. /* given the mean anomaly, ma, and the eccentricity, s, of elliptical motion,
  154.  * find the true anomaly, *nu, and the eccentric anomaly, *ea.
  155.  * all angles in radians.
  156.  */
  157. anomaly (ma, s, nu, ea)
  158. float ma, s;
  159. float *nu, *ea;
  160. {
  161.     float m, dla, fea;
  162.  
  163.     m = ma-TWOPI*(int)(ma/TWOPI);
  164.     fea = m;
  165.     while (1) {
  166.         dla = fea-(s*sin(fea))-m;
  167.         if (fabs(dla)<1e-6)
  168.         break;
  169.         dla /= 1-(s*cos(fea));
  170.         fea -= dla;
  171.     }
  172.  
  173.     *nu = 2*atan(sqrt((1+s)/(1-s))*tan(fea/2));
  174.     *ea = fea;
  175. }
  176. xXx
  177. echo extracting astro.h
  178. cat > astro.h << 'xXx'
  179. #ifndef PI
  180. #define    PI        3.141592653589793
  181. #endif
  182.  
  183. /* conversions among hours (of ra), degrees and radians. */
  184. #define    degrad(x)    ((x)*PI/180.)
  185. #define    raddeg(x)    ((x)*180./PI)
  186. #define    hrdeg(x)    ((x)*15.)
  187. #define    deghr(x)    ((x)/15.)
  188. #define    hrrad(x)    degrad(hrdeg(x))
  189. #define    radhr(x)    deghr(raddeg(x))
  190.  
  191. /* manifest names for planets.
  192.  * N.B. must cooincide with usage in pelement.c and plans.c.
  193.  */
  194. #define    MERCURY    0
  195. #define    VENUS    1
  196. #define    MARS    2
  197. #define    JUPITER    3
  198. #define    SATURN    4
  199. #define    URANUS    5
  200. #define    NEPTUNE    6
  201. #define    PLUTO    7
  202. xXx
  203. echo extracting cal_mjd.c
  204. cat > cal_mjd.c << 'xXx'
  205. #include <stdio.h>
  206. #include <math.h>
  207. #include "astro.h"
  208.  
  209. /* given a date in months, mn, days, dy, years, yr,
  210.  * return the modified Julian date (number of days elapsed since 1900 jan 0.5),
  211.  * *mjd.
  212.  */
  213. cal_mjd (mn, dy, yr, mjd)
  214. int mn, yr;
  215. float dy;
  216. float *mjd;
  217. {
  218.     int b, d, m, y;
  219.     long c;
  220.  
  221.     m = mn;
  222.     y = (yr < 0) ? yr + 1 : yr;
  223.     if (mn < 3) {
  224.         m += 12;
  225.         y -= 1;
  226.     }
  227.  
  228.     if (yr < 1582 || yr == 1582 && (mn < 10 || mn == 10 && dy < 15)) 
  229.         b = 0;
  230.     else {
  231.         int a;
  232.         a = y/100;
  233.         b = 2 - a + a/4;
  234.     }
  235.  
  236.     if (y < 0)
  237.         c = (long)((365.25*y) - 0.75) - 694025L;
  238.     else
  239.         c = (long)(365.25*y) - 694025L;
  240.  
  241.     d = 30.6001*(m+1);
  242.  
  243.     *mjd = b + c + d + dy - 0.5;
  244. }
  245.  
  246. /* given the modified Julian date (number of days elapsed since 1900 jan 0.5,),
  247.  * mjd, return the calendar date in months, *mn, days, *dy, and years, *yr.
  248.  */
  249. mjd_cal (mjd, mn, dy, yr)
  250. float mjd;
  251. int *mn, *yr;
  252. float *dy;
  253. {
  254.     float d, f;
  255.     float i, a, b, ce, g;
  256.  
  257.     d = mjd + 0.5;
  258.     i = floor(d);
  259.     f = d-i;
  260.     if (f == 1) {
  261.         f = 0;
  262.         i += 1;
  263.     }
  264.  
  265.     if (i > -115860.0) {
  266.         a = floor((i/36524.25)+.9983573)+14;
  267.         i += 1 + a - floor(a/4.0);
  268.     }
  269.  
  270.     b = floor((i/365.25)+.802601);
  271.     ce = i - floor((365.25*b)+.750001)+416;
  272.     g = floor(ce/30.6001);
  273.     *mn = g - 1;
  274.     *dy = ce - floor(30.6001*g)+f;
  275.     *yr = b + 1899;
  276.  
  277.     if (g > 13.5)
  278.         *mn = g - 13;
  279.     if (*mn < 2.5)
  280.         *yr = b + 1900;
  281.     if (*yr < 1)
  282.         *yr -= 1;
  283. }
  284.  
  285. /* given an mjd, set *dow to 0..6 according to which dayof the week it falls
  286.  * on (0=sunday) or set it to -1 if can't figure it out.
  287.  */
  288. mjd_dow (mjd, dow)
  289. float mjd;
  290. int *dow;
  291. {
  292.     /* cal_mjd() uses Gregorian dates on or after Oct 15, 1582.
  293.      * (Pope Gregory XIII dropped 10 days, Oct 5..14, and improved the leap-
  294.      * year algorithm). however, Great Britian and the colonies did not
  295.      * adopt it until Sept 14, 1752 (they dropped 11 days, Sept 3-13,
  296.      * due to additional accumulated error). leap years before 1752 thus
  297.      * can not easily be accounted for from the cal_mjd() number...
  298.      */
  299.     if (mjd < -53798.5) {
  300.         /* pre sept 14, 1752 too hard to correct */
  301.         *dow = -1;
  302.         return;
  303.     }
  304.     *dow = ((int)floor(mjd-.5) + 1) % 7; /* 1/1/1900 (mjd 0.5) is a Monday*/
  305.     if (*dow < 0)
  306.         *dow += 7;
  307. }
  308.  
  309. /* given a mjd, return the the number of days in the month.  */
  310. mjd_dpm (mjd, ndays)
  311. float mjd;
  312. int *ndays;
  313. {
  314.     static short dpm[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
  315.     int m, y;
  316.     float d;
  317.  
  318.     mjd_cal (mjd, &m, &d, &y);
  319.     *ndays = (m==2 && ((y%4==0 && y%100!=0)||y%400==0)) ? 29 : dpm[m-1];
  320. }
  321. xXx
  322. echo extracting eq_ecl.c 
  323. cat > eq_ecl.c << 'xXx'
  324. #include <stdio.h>
  325. #include <math.h>
  326. #include "astro.h"
  327.  
  328. #define    EQtoECL    1
  329. #define    ECLtoEQ    (-1)
  330.  
  331. /* given the modified Julian date, mjd, and an equitorial ra and dec, each in
  332.  * radians, find the corresponding geocentric ecliptic latitude, *lat, and
  333.  * longititude, *lng, also each in radians.
  334.  * the effects of nutation and the angle of the obliquity are included.
  335.  */
  336. eq_ecl (mjd, ra, dec, lat, lng)
  337. float mjd, ra, dec;
  338. float *lat, *lng;
  339. {
  340.     ecleq_aux (EQtoECL, mjd, ra, dec, lng, lat);
  341. }
  342.  
  343. /* given the modified Julian date, mjd, and a geocentric ecliptic latitude,
  344.  * *lat, and longititude, *lng, each in radians, find the corresponding
  345.  * equitorial ra and dec, also each in radians.
  346.  * the effects of nutation and the angle of the obliquity are included.
  347.  */
  348. ecl_eq (mjd, lat, lng, ra, dec)
  349. float mjd, lat, lng;
  350. float *ra, *dec;
  351. {
  352.     ecleq_aux (ECLtoEQ, mjd, lng, lat, ra, dec);
  353. }
  354.  
  355. static
  356. ecleq_aux (sw, mjd, x, y, p, q)
  357. int sw;            /* +1 for eq to ecliptic, -1 for vv. */
  358. float mjd, x, y;    /* sw==1: x==ra, y==dec.  sw==-1: x==lng, y==lat. */
  359. float *p, *q;        /* sw==1: p==lng, q==lat. sw==-1: p==ra, q==dec. */
  360. {
  361.     static float lastmjd;        /* last mjd calculated */
  362.     static float seps, ceps;    /* sin and cos of mean obliquity */
  363.     float sx, cx, sy, cy, ty;
  364.  
  365.     if (mjd != lastmjd) {
  366.         float deps, dpsi, eps;
  367.         obliquity (mjd, &eps);        /* mean obliquity for date */
  368. #ifdef NONUTATION
  369. #define NONUTATION
  370.         deps = 0;                /* checkout handler did not
  371.                          * include nutation correction.
  372.                          */
  373. #else
  374.         nutation (mjd, &deps, &dpsi);    /* correctin due to nutation */
  375. #endif
  376.         eps += deps;            /* true obliquity for date */
  377.             seps = sin(eps);
  378.         ceps = cos(eps);
  379.         lastmjd = mjd;
  380.     }
  381.  
  382.     sy = sin(y);
  383.     cy = cos(y);                /* always non-negative */
  384.         if (fabs(cy)<1e-20) cy = 1e-20;        /* insure > 0 */
  385.         ty = sy/cy;
  386.     cx = cos(x);
  387.     sx = sin(x);
  388.         *q = asin((sy*ceps)-(cy*seps*sx*sw));
  389.         *p = atan(((sx*ceps)+(ty*seps*sw))/cx);
  390.         if (cx<0) *p += PI;        /* account for atan quad ambiguity */
  391.     range (p, 2*PI);
  392. }
  393. xXx
  394. echo extracting moon.c 
  395. cat > moon.c << 'xXx'
  396. #include <stdio.h>
  397. #include <math.h>
  398. #include "astro.h"
  399.  
  400. /* given the mjd, find the geocentric ecliptic longitude, lam, and latitude,
  401.  * bet, and horizontal parallax, hp for the moon.
  402.  * N.B. series for long and lat are good to about 10 and 3 arcseconds. however,
  403.  *   math errors cause up to 100 and 30 arcseconds error, even if use double.
  404.  *   why?? suspect highly sensitive nature of difference used to get m1..6.
  405.  * N.B. still need to correct for nutation. then for topocentric location
  406.  *   further correct for parallax and refraction.
  407.  */
  408. moon (mjd, lam, bet, hp)
  409. float mjd;
  410. float *lam, *bet, *hp;
  411. {
  412.     float t, t2;
  413.     float ld;
  414.     float ms;
  415.     float md;
  416.     float de;
  417.     float f;
  418.     float n;
  419.     float a, sa, sn, b, sb, c, sc, e, e2, l, g, w1, w2;
  420.     float m1, m2, m3, m4, m5, m6;
  421.  
  422.     t = mjd/36525.;
  423.     t2 = t*t;
  424.  
  425.     m1 = mjd/27.32158213;
  426.     m1 = 360.0*(m1-(int)m1);
  427.     m2 = mjd/365.2596407;
  428.     m2 = 360.0*(m2-(int)m2);
  429.     m3 = mjd/27.55455094;
  430.     m3 = 360.0*(m3-(int)m3);
  431.     m4 = mjd/29.53058868;
  432.     m4 = 360.0*(m4-(int)m4);
  433.     m5 = mjd/27.21222039;
  434.     m5 = 360.0*(m5-(int)m5);
  435.     m6 = mjd/6798.363307;
  436.     m6 = 360.0*(m6-(int)m6);
  437.  
  438.     ld = 270.434164+m1-(.001133-.0000019*t)*t2;
  439.     ms = 358.475833+m2-(.00015+.0000033*t)*t2;
  440.     md = 296.104608+m3+(.009192+.0000144*t)*t2;
  441.     de = 350.737486+m4-(.001436-.0000019*t)*t2;
  442.     f = 11.250889+m5-(.003211+.0000003*t)*t2;
  443.     n = 259.183275-m6+(.002078+.000022*t)*t2;
  444.  
  445.     a = degrad(51.2+20.2*t);
  446.     sa = sin(a);
  447.     sn = sin(degrad(n));
  448.     b = 346.56+(132.87-.0091731*t)*t;
  449.     sb = .003964*sin(degrad(b));
  450.     c = degrad(n+275.05-2.3*t);
  451.     sc = sin(c);
  452.     ld = ld+.000233*sa+sb+.001964*sn;
  453.     ms = ms-.001778*sa;
  454.     md = md+.000817*sa+sb+.002541*sn;
  455.     f = f+sb-.024691*sn-.004328*sc;
  456.     de = de+.002011*sa+sb+.001964*sn;
  457.     e = 1-(.002495+7.52e-06*t)*t;
  458.     e2 = e*e;
  459.  
  460.     ld = degrad(ld);
  461.     ms = degrad(ms);
  462.     n = degrad(n);
  463.     de = degrad(de);
  464.     f = degrad(f);
  465.     md = degrad(md);
  466.  
  467.     l = 6.28875*sin(md)+1.27402*sin(2*de-md)+.658309*sin(2*de)+
  468.         .213616*sin(2*md)-e*.185596*sin(ms)-.114336*sin(2*f)+
  469.         .058793*sin(2*(de-md))+.057212*e*sin(2*de-ms-md)+
  470.         .05332*sin(2*de+md)+.045874*e*sin(2*de-ms)+.041024*e*sin(md-ms);
  471.     l = l-.034718*sin(de)-e*.030465*sin(ms+md)+.015326*sin(2*(de-f))-
  472.         .012528*sin(2*f+md)-.01098*sin(2*f-md)+.010674*sin(4*de-md)+
  473.         .010034*sin(3*md)+.008548*sin(4*de-2*md)-e*.00791*sin(ms-md+2*de)-
  474.         e*.006783*sin(2*de+ms);
  475.     l = l+.005162*sin(md-de)+e*.005*sin(ms+de)+.003862*sin(4*de)+
  476.         e*.004049*sin(md-ms+2*de)+.003996*sin(2*(md+de))+
  477.         .003665*sin(2*de-3*md)+e*.002695*sin(2*md-ms)+
  478.         .002602*sin(md-2*(f+de))+e*.002396*sin(2*(de-md)-ms)-
  479.         .002349*sin(md+de);
  480.     l = l+e2*.002249*sin(2*(de-ms))-e*.002125*sin(2*md+ms)-
  481.         e2*.002079*sin(2*ms)+e2*.002059*sin(2*(de-ms)-md)-
  482.         .001773*sin(md+2*(de-f))-.001595*sin(2*(f+de))+
  483.         e*.00122*sin(4*de-ms-md)-.00111*sin(2*(md+f))+.000892*sin(md-3*de);
  484.     l = l-e*.000811*sin(ms+md+2*de)+e*.000761*sin(4*de-ms-2*md)+
  485.          e2*.000704*sin(md-2*(ms+de))+e*.000693*sin(ms-2*(md-de))+
  486.          e*.000598*sin(2*(de-f)-ms)+.00055*sin(md+4*de)+.000538*sin(4*md)+
  487.          e*.000521*sin(4*de-ms)+.000486*sin(2*md-de);
  488.     l = l+e2*.000717*sin(md-2*ms);
  489.     *lam = ld+degrad(l);
  490.     range (lam, 2*PI);
  491.  
  492.     g = 5.12819*sin(f)+.280606*sin(md+f)+.277693*sin(md-f)+
  493.         .173238*sin(2*de-f)+.055413*sin(2*de+f-md)+.046272*sin(2*de-f-md)+
  494.         .032573*sin(2*de+f)+.017198*sin(2*md+f)+.009267*sin(2*de+md-f)+
  495.         .008823*sin(2*md-f)+e*.008247*sin(2*de-ms-f);
  496.     g = g+.004323*sin(2*(de-md)-f)+.0042*sin(2*de+f+md)+
  497.         e*.003372*sin(f-ms-2*de)+e*.002472*sin(2*de+f-ms-md)+
  498.         e*.002222*sin(2*de+f-ms)+e*.002072*sin(2*de-f-ms-md)+
  499.         e*.001877*sin(f-ms+md)+.001828*sin(4*de-f-md)-e*.001803*sin(f+ms)-
  500.         .00175*sin(3*f);
  501.     g = g+e*.00157*sin(md-ms-f)-.001487*sin(f+de)-e*.001481*sin(f+ms+md)+
  502.          e*.001417*sin(f-ms-md)+e*.00135*sin(f-ms)+.00133*sin(f-de)+
  503.          .001106*sin(f+3*md)+.00102*sin(4*de-f)+.000833*sin(f+4*de-md)+
  504.          .000781*sin(md-3*f)+.00067*sin(f+4*de-2*md);
  505.     g = g+.000606*sin(2*de-3*f)+.000597*sin(2*(de+md)-f)+
  506.         e*.000492*sin(2*de+md-ms-f)+.00045*sin(2*(md-de)-f)+
  507.         .000439*sin(3*md-f)+.000423*sin(f+2*(de+md))+
  508.         .000422*sin(2*de-f-3*md)-e*.000367*sin(ms+f+2*de-md)-
  509.         e*.000353*sin(ms+f+2*de)+.000331*sin(f+4*de);
  510.     g = g+e*.000317*sin(2*de+f-ms+md)+e2*.000306*sin(2*(de-ms)-f)-
  511.         .000283*sin(md+3*f);
  512.     w1 = .0004664*cos(n);
  513.     w2 = .0000754*cos(c);
  514.     *bet = degrad(g)*(1-w1-w2);
  515.  
  516.     *hp = .950724+.051818*cos(md)+.009531*cos(2*de-md)+.007843*cos(2*de)+
  517.           .002824*cos(2*md)+.000857*cos(2*de+md)+e*.000533*cos(2*de-ms)+
  518.           e*.000401*cos(2*de-md-ms)+e*.00032*cos(md-ms)-.000271*cos(de)-
  519.           e*.000264*cos(ms+md)-.000198*cos(2*f-md);
  520.     *hp = *hp+.000173*cos(3*md)+.000167*cos(4*de-md)-e*.000111*cos(ms)+
  521.          .000103*cos(4*de-2*md)-.000084*cos(2*md-2*de)-
  522.          e*.000083*cos(2*de+ms)+.000079*cos(2*de+2*md)+.000072*cos(4*de)+
  523.          e*.000064*cos(2*de-ms+md)-e*.000063*cos(2*de+ms-md)+
  524.          e*.000041*cos(ms+de);
  525.     *hp = *hp+e*.000035*cos(2*md-ms)-.000033*cos(3*md-2*de)-
  526.          .00003*cos(md+de)-.000029*cos(2*(f-de))-e*.000029*cos(2*md+ms)+
  527.          e2*.000026*cos(2*(de-ms))-.000023*cos(2*(f-de)+md)+
  528.          e*.000019*cos(4*de-ms-md);
  529.     *hp = degrad(*hp);
  530. }
  531. xXx
  532. echo extracting moonnf.c
  533. cat > moonnf.c << 'xXx'
  534. #include <stdio.h>
  535. #include <math.h>
  536. #include "astro.h"
  537.  
  538. #define    unw(w,z)    ((w)-floor((w)/(z))*(z))
  539.  
  540. /* given a modified Julian date, mjd, return the mjd of the new
  541.  * and full moons about then, mjdn and mjdf.
  542.  * TODO: exactly which ones does it find? eg:
  543.  *   5/28/1988 yields 5/15 and 5/31
  544.  *   5/29             6/14     6/29
  545.  */
  546. moonnf (mjd, mjdn, mjdf)
  547. float mjd;
  548. float *mjdn, *mjdf;
  549. {
  550.     int mo, yr;
  551.     float dy;
  552.     float mjd0;
  553.     float k, tn, tf, t;
  554.  
  555.     mjd_cal (mjd, &mo, &dy, &yr);
  556.     cal_mjd (1, 0., yr, &mjd0);
  557.     k = (yr-1900+((mjd-mjd0)/365))*12.3685;
  558.     k = floor(k+0.5);
  559.     tn = k/1236.85;
  560.     tf = (k+0.5)/1236.85;
  561.     t = tn;
  562.     m (t, k, mjdn);
  563.     t = tf;
  564.     k += 0.5;
  565.     m (t, k, mjdf);
  566. }
  567.  
  568. static
  569. m (t, k, mjd)
  570. float t, k;
  571. float *mjd;
  572. {
  573.     float t2, a, a1, b, b1, c, ms, mm, f, ddjd;
  574.  
  575.     t2 = t*t;
  576.     a = 29.53*k;
  577.     c = degrad(166.56+(132.87-9.173e-3*t)*t);
  578.     b = 5.8868e-4*k+(1.178e-4-1.55e-7*t)*t2+3.3e-4*sin(c)+7.5933E-1;
  579.     ms = 359.2242+360*unw(k/1.236886e1,1)-(3.33e-5+3.47e-6*t)*t2;
  580.     mm = 306.0253+360*unw(k/9.330851e-1,1)+(1.07306e-2+1.236e-5*t)*t2;
  581.     f = 21.2964+360*unw(k/9.214926e-1,1)-(1.6528e-3+2.39e-6*t)*t2;
  582.     ms = unw(ms,360);
  583.     mm = unw(mm,360);
  584.     f = unw(f,360);
  585.     ms = degrad(ms);
  586.     mm = degrad(mm);
  587.     f = degrad(f);
  588.     ddjd = (1.734e-1-3.93e-4*t)*sin(ms)+2.1e-3*sin(2*ms)
  589.         -4.068e-1*sin(mm)+1.61e-2*sin(2*mm)-4e-4*sin(3*mm)
  590.         +1.04e-2*sin(2*f)-5.1e-3*sin(ms+mm)-7.4e-3*sin(ms-mm)
  591.         +4e-4*sin(2*f+ms)-4e-4*sin(2*f-ms)-6e-4*sin(2*f+mm)
  592.         +1e-3*sin(2*f-mm)+5e-4*sin(ms+2*mm);
  593.     a1 = (int)a;
  594.     b = b+ddjd+(a-a1);
  595.     b1 = (int)b;
  596.     a = a1+b1;
  597.     b = b-b1;
  598.     *mjd = a + b;
  599. }
  600. xXx
  601. echo extracting moonrs.c
  602. cat > moonrs.c << 'xXx'
  603. #include <stdio.h>
  604. #include <math.h>
  605. #include "astro.h"
  606.  
  607. #define    PASSES    2
  608.  
  609. /* given the mjd and geographical latitude and longitude, lat and lng, 
  610.  *   find the ut and azimuth of moonrise, utcr and azr, and moonset, utcs and
  611.  *   azs for the day of mjd. all angles are in radians.
  612.  * times are for the upper limb, including effects of refraction (a constant 34
  613.  *   arcminutes), parallax, and nutation. accuracy is to nearest minute.
  614.  * the calculated times can be in error if they occur within about a half hour 
  615.  *   of utc midnight due to the algorithm's confusion about the date. the best
  616.  *   way to verify it is to calculate the times for several days either side of
  617.  *   the date to see whether the problem values conform to a smooth trend set
  618.  *   by the others.
  619.  * possible values of status:
  620.  *    2: can't cope (such as geographical latitude very near +-90).
  621.  *    1: moon never rises/sets on date.
  622.  *    0: all ok.
  623.  *   -1: moon is circumpolar.
  624.  *   -2: possible error in near-midnight moonrise time; see above.
  625.  *   -3: possible error in near-midnight moonset time; see above.
  626.  */
  627. moonrs (mjd, lat, lng, utcr, utcs, azr, azs, status)
  628. float mjd, lat, lng;
  629. float *utcr, *utcs;
  630. float *azr, *azs;
  631. int *status;
  632. {
  633.     float lstr, lsts;    /* local sidereal times of rising/setting */
  634.     float al = radhr(lng);    /* time correction for longitude */
  635.     int midnight = 0;    /* midnight caution */
  636.     float mjd0;        /* start of mjd day */
  637.     float x;        /* discarded tmp value */
  638.     float ut;
  639.     int pass;
  640.  
  641.     mjd0 = floor(mjd+0.5)-0.5;
  642.  
  643.     /* first approximation is to find rise/set times of a fixed object
  644.      * in the position of the moon at local midday.
  645.      */
  646.     fixedmoonriset (mjd0+(12.0-al)/24.0, lat, &lstr, &lsts, &x, &x, status);
  647.     if (*status != 0) return;
  648.  
  649.     /* find a better approximation to the rising circumstances based on two
  650.      * more passes, each using a "fixed" object at the moons location at
  651.      * previous approximation of the rise time.
  652.      */
  653.     pass = 0;
  654.     while (1) {
  655.         gst_utc (mjd0, lstr, &ut);
  656.         if (++pass > PASSES)
  657.         break;
  658.         fixedmoonriset (mjd0+(ut-al)/24., lat, &lstr, &x, azr, &x, status);
  659.         if (*status != 0) return;
  660.     }
  661.     if (ut > 23.5 || ut < 0.5)
  662.         midnight = -2;        /* moonrise caution near midnight */
  663.     *utcr = ut - (al * .9972696);    /* allow for sidereal shift */
  664.  
  665.     /* find a better approximation to the setting circumstances based on two
  666.      * more passes, each using a "fixed" object at the moons location at
  667.      * previous approximation of the setting time.
  668.      */
  669.     pass = 0;
  670.     while (1) {
  671.         gst_utc (mjd0, lsts, &ut);
  672.         if (++pass > PASSES)
  673.         break;
  674.         fixedmoonriset (mjd0+(ut-al)/24., lat, &x, &lsts, &x, azs, status);
  675.         if (*status != 0) return;
  676.     }
  677.     if (ut > 23.5 || ut < 0.5)
  678.         midnight = -3;        /* moonset caution near midnight */
  679.     *utcs = ut - (al * .9972696);    /* allow for sidereal shift */
  680.  
  681.     if (midnight)
  682.         *status = midnight;        /* report caution if near midnight */
  683. }
  684.  
  685. /* find the local rise/set sidereal times and azimuths of an object fixed
  686.  * at the ra/dec of the moon on the given mjd time as seen from lat.
  687.  */
  688. static
  689. fixedmoonriset (mjd, lat, lstr, lsts, azr, azs, status)
  690. float mjd, lat;
  691. float *lstr, *lsts;
  692. float *azr, *azs;
  693. int *status;
  694. {
  695.     float lam, bet, hp;
  696.     float deps, dpsi;
  697.     float dis;
  698.     float ra, dec;
  699.  
  700.     /* find geocentric ecliptic location and equatorial horizontal parallax
  701.      * of moon at mjd.
  702.      */
  703.     moon (mjd, &lam, &bet, &hp);
  704.  
  705.     /* allow for nutation */
  706.     nutation (mjd, &deps, &dpsi);
  707.         lam += dpsi;
  708.  
  709.     /* convert ecliptic to equatorial coords */
  710.     ecl_eq (mjd, bet, lam, &ra, &dec);
  711.  
  712.     /* find local sidereal times of rise/set times/azimuths for given
  713.      * equatorial coords, allowing for
  714.      * .27249*sin(hp)   parallax (hp is radius of earth from moon;
  715.      *                  .27249 is radius of moon from earth where
  716.      *                  the ratio is the dia_moon/dia_earth).
  717.      * 34' (9.8902e-3)  nominal refraction
  718.      * hp               nominal angular earth radius. subtract because
  719.      *                  tangential line-of-sight makes moon appear lower
  720.      *                  as move out from center of earth.
  721.      * TODO: do better at refraction; see fixedsunriset() in sunrs.c.
  722.      */
  723.     dis = .27249*sin(hp) + 9.8902e-3 - hp;
  724.     riset (ra, dec, lat, dis, lstr, lsts, azr, azs, status);
  725. }
  726. xXx
  727. echo extracting nutation.c
  728. cat > nutation.c << 'xXx'
  729. #include <stdio.h>
  730. #include <math.h>
  731. #include "astro.h"
  732.  
  733. /* given the modified JD, mjd, find the nutation in obliquity, *deps, and
  734.  * the nutation in longitude, *dpsi, each in radians.
  735.  */
  736. nutation (mjd, deps, dpsi)
  737. float mjd;
  738. float *deps, *dpsi;
  739. {
  740.     static float lastmjd, lastdeps, lastdpsi;
  741.     float ls, ld;    /* sun's mean longitude, moon's mean longitude */
  742.     float ms, md;    /* sun's mean anomaly, moon's mean anomaly */
  743.     float nm;    /* longitude of moon's ascending node */
  744.     float t, t2;    /* number of Julian centuries of 36525 days since
  745.              * Jan 0.5 1900.
  746.              */
  747.     float tls, tnm, tld;    /* twice above */
  748.     double a, b;    /* temps */
  749.  
  750.     if (mjd == lastmjd) {
  751.         *deps = lastdeps;
  752.         *dpsi = lastdpsi;
  753.         return;
  754.     }
  755.         
  756.     t = mjd/36525.;
  757.     t2 = t*t;
  758.  
  759.     a = 100.0021358*t;
  760.     b = 360.*(a-(int)a);
  761.     ls = 279.697+.000303*t2+b;
  762.  
  763.     a = 1336.855231*t;
  764.     b = 360.*(a-(int)a);
  765.     ld = 270.434-.001133*t2+b;
  766.  
  767.     a = 99.99736056000026*t;
  768.     b = 360.*(a-(int)a);
  769.     ms = 358.476-.00015*t2+b;
  770.  
  771.     a = 13255523.59*t;
  772.     b = 360.*(a-(int)a);
  773.     md = 296.105+.009192*t2+b;
  774.  
  775.     a = 5.372616667*t;
  776.     b = 360.*(a-(int)a);
  777.     nm = 259.183+.002078*t2-b;
  778.  
  779.     /* convert to radian forms for use with trig functions.
  780.      */
  781.     tls = 2*degrad(ls);
  782.     nm = degrad(nm);
  783.     tnm = 2*degrad(nm);
  784.     ms = degrad(ms);
  785.     tld = 2*degrad(ld);
  786.     md = degrad(md);
  787.  
  788.     /* find delta psi and eps, in arcseconds.
  789.      */
  790.     lastdpsi = (-17.2327-.01737*t)*sin(nm)+(-1.2729-.00013*t)*sin(tls)
  791.            +.2088*sin(tnm)-.2037*sin(tld)+(.1261-.00031*t)*sin(ms)
  792.            +.0675*sin(md)-(.0497-.00012*t)*sin(tls+ms)
  793.            -.0342*sin(tld-nm)-.0261*sin(tld+md)+.0214*sin(tls-ms)
  794.            -.0149*sin(tls-tld+md)+.0124*sin(tls-nm)+.0114*sin(tld-md);
  795.     lastdeps = (9.21+.00091*t)*cos(nm)+(.5522-.00029*t)*cos(tls)
  796.            -.0904*cos(tnm)+.0884*cos(tld)+.0216*cos(tls+ms)
  797.            +.0183*cos(tld-nm)+.0113*cos(tld+md)-.0093*cos(tls-ms)
  798.            -.0066*cos(tls-nm);
  799.  
  800.     /* convert to radians.
  801.      */
  802.     lastdpsi = degrad(lastdpsi/3600);
  803.     lastdeps = degrad(lastdeps/3600);
  804.  
  805.     lastmjd = mjd;
  806.     *deps = lastdeps;
  807.     *dpsi = lastdpsi;
  808. }
  809. xXx
  810. echo extracting obliq.c
  811. cat > obliq.c << 'xXx'
  812. #include <stdio.h>
  813. #include "astro.h"
  814.  
  815. /* given the modified Julian date, mjd, find the obliquity of the
  816.  * ecliptic, *eps, in radians.
  817.  */
  818. obliquity (mjd, eps)
  819. float mjd;
  820. float *eps;
  821. {
  822.     static float lastmjd, lasteps;
  823.  
  824.     if (mjd != lastmjd) {
  825.         float t;
  826.         t = mjd/36525.;
  827.         lasteps = degrad(2.345229444E1 - ((((-1.81E-3*t)+5.9E-3)*t+4.6845E1)*t)/3600.0);
  828.         lastmjd = mjd;
  829.     }
  830.     *eps = lasteps;
  831. }
  832. xXx
  833. echo extracting parallax.c
  834. cat > parallax.c << 'xXx'
  835. #include <stdio.h>
  836. #include <math.h>
  837. #include "astro.h"
  838.  
  839. /* given true ha and dec, tha and tdec, the geographical latitude, phi, the
  840.  * height above sea-level (as a fraction of the earths radius, 6378.16km),
  841.  * ht, and the equatorial horizontal parallax, ehp, find the apparent
  842.  * ha and dec, aha and adec allowing for parallax.
  843.  * all angles in radians. ehp is the angle subtended at the body by the
  844.  * earth's equator.
  845.  */
  846. ta_par (tha, tdec, phi, ht, ehp, aha, adec)
  847. float tha, tdec, phi, ht, ehp;
  848. float *aha, *adec;
  849. {
  850.     static float last_phi, last_ht, rsp, rcp;
  851.     float rp;    /* distance to object in Earth radii */
  852.     float ctha;
  853.     float stdec, ctdec;
  854.     float tdtha, dtha;
  855.     float caha;
  856.  
  857.     /* avoid calcs involving the same phi and ht */
  858.     if (phi != last_phi || ht != last_ht) {
  859.         float cphi, sphi, u;
  860.         cphi = cos(phi);
  861.         sphi = sin(phi);
  862.         u = atan(9.96647e-1*sphi/cphi);
  863.         rsp = (9.96647e-1*sin(u))+(ht*sphi);
  864.         rcp = cos(u)+(ht*cphi);
  865.         last_phi  =  phi;
  866.         last_ht  =  ht;
  867.     }
  868.  
  869.         rp = 1/sin(ehp);
  870.  
  871.         ctha = cos(tha);
  872.     stdec = sin(tdec);
  873.     ctdec = cos(tdec);
  874.         tdtha = (rcp*sin(tha))/((rp*ctdec)-(rcp*ctha));
  875.         dtha = atan(tdtha);
  876.     *aha = tha+dtha;
  877.     caha = cos(*aha);
  878.     range (aha, 2*PI);
  879.         *adec = atan(caha*(rp*stdec-rsp)/(rp*ctdec*ctha-rcp));
  880. }
  881.  
  882. /* given the apparent ha and dec, aha and adec, the geographical latitude, phi,
  883.  * the height above sea-level (as a fraction of the earths radius, 6378.16km),
  884.  * ht, and the equatorial horizontal parallax, ehp, find the true ha and dec,
  885.  * tha and tdec allowing for parallax.
  886.  * all angles in radians. ehp is the angle subtended at the body by the
  887.  * earth's equator.
  888.  * uses ta_par() iteratively: find a set of true ha/dec that converts back
  889.   *  to the given apparent ha/dec.
  890.  */
  891. at_par (aha, adec, phi, ht, ehp, tha, tdec)
  892. float aha, adec, phi, ht, ehp;
  893. float *tha, *tdec;
  894. {
  895.     float nha, ndec;    /* ha/dec corres. to current true guesses */
  896.     float eha, edec;    /* error in ha/dec */
  897.  
  898.     /* first guess for true is just the apparent */
  899.     *tha = aha;
  900.     *tdec = adec;
  901.  
  902.     while (1) {
  903.         ta_par (*tha, *tdec, phi, ht, ehp, &nha, &ndec);
  904.         eha = aha - nha;
  905.         edec = adec - ndec;
  906.         if (fabs(eha)<1e-6 && fabs(edec)<1e-6)
  907.         break;
  908.         *tha += eha;
  909.         *tdec += edec;
  910.     }
  911. }
  912. xXx
  913. echo extracting pelement.c
  914. cat > pelement.c << 'xXx'
  915. #include <stdio.h>
  916. #include <math.h>
  917. #include "astro.h"
  918.  
  919. /* this array contains polynomial coefficients to find the various orbital
  920.  *   elements for the mean orbit at any instant in time for each major planet.
  921.  *   the first five elements are in the form a0 + a1*t + a2*t**2 + a3*t**3,
  922.  *   where t is the number of Julian centuries of 36525 Julian days since 1900
  923.  *   Jan 0.5. the last three elements are constants.
  924.  *
  925.  * the orbital element (column) indeces are:
  926.  *   [ 0- 3]: coefficients for mean longitude, in degrees;
  927.  *   [ 4- 7]: coefficients for longitude of the perihelion, in degrees;
  928.  *   [ 8-11]: coefficients for eccentricity;
  929.  *   [12-15]: coefficients for inclination, in degrees;
  930.  *   [16-19]: coefficients for longitude of the ascending node, in degrees;
  931.  *      [20]: semi-major axis, in AU;
  932.  *      [21]: angular diameter at 1 AU, in arcsec;
  933.  *      [22]: standard visual magnitude, ie, the visual magnitude of the planet
  934.  *          when at a distance of 1 AU from both the Sun and the Earth and
  935.  *          with zero phase angle.
  936.  *
  937.  * the planent (row) indeces are:
  938.  *   [0]: Mercury; [1]: Venus;   [2]: Mars;  [3]: Jupiter; [4]: Saturn;
  939.  *   [5]: Uranus;  [6]: Neptune; [7]: Pluto.
  940.  */
  941. #define    NPELE    (5*4 + 3)    /* 4 coeffs for ea of 5 elems, + 3 constants */
  942. static float elements[8][NPELE] = {
  943.  
  944.     {   /*     mercury... */
  945.  
  946.         178.179078,    415.2057519,    3.011e-4,    0.0,
  947.         75.899697,    1.5554889,    2.947e-4,    0.0,
  948.         .20561421,    2.046e-5,    3e-8,        0.0,
  949.         7.002881,    1.8608e-3,    -1.83e-5,    0.0,
  950.         47.145944,    1.1852083,    1.739e-4,    0.0,
  951.         .3870986,    6.74,         -0.42
  952.     },
  953.  
  954.     {   /*     venus... */
  955.  
  956.         342.767053,    162.5533664,    3.097e-4,    0.0,
  957.         130.163833,    1.4080361,    -9.764e-4,    0.0,
  958.         6.82069e-3,    -4.774e-5,    9.1e-8,        0.0,
  959.         3.393631,    1.0058e-3,    -1e-6,        0.0,
  960.         75.779647,    .89985,        4.1e-4,        0.0,
  961.         .7233316,    16.92,        -4.4
  962.     },
  963.  
  964.     {   /*     mars... */
  965.  
  966.         293.737334,    53.17137642,    3.107e-4,    0.0,
  967.         3.34218203e2, 1.8407584,    1.299e-4,    -1.19e-6,
  968.         9.33129e-2,    9.2064e-5,    7.7e-8,        0.0,
  969.         1.850333,    -6.75e-4,    1.26e-5,    0.0,
  970.         48.786442,    .7709917,    -1.4e-6,    -5.33e-6,
  971.         1.5236883,    9.36,        -1.52
  972.     },
  973.  
  974.     {   /*     jupiter... */
  975.  
  976.         238.049257,    8.434172183,    3.347e-4,    -1.65e-6,
  977.         1.2720972e1, 1.6099617,    1.05627e-3,    -3.43e-6,
  978.         4.833475e-2, 1.6418e-4,    -4.676e-7,    -1.7e-9,
  979.         1.308736,    -5.6961e-3,    3.9e-6,        0.0,
  980.         99.443414,    1.01053,    3.5222e-4,    -8.51e-6,
  981.         5.202561,    196.74,        -9.4
  982.     },
  983.  
  984.     {   /*     saturn... */
  985.  
  986.         266.564377,    3.398638567,    3.245e-4,    -5.8e-6,
  987.         9.1098214e1, 1.9584158,    8.2636e-4,    4.61e-6,
  988.         5.589232e-2, -3.455e-4,    -7.28e-7,    7.4e-10,
  989.         2.492519,    -3.9189e-3,    -1.549e-5,    4e-8,
  990.         112.790414,    .8731951,    -1.5218e-4,    -5.31e-6,
  991.         9.554747,    165.6,        -8.88
  992.     },
  993.  
  994.     {   /*     uranus... */
  995.  
  996.         244.19747,    1.194065406,    3.16e-4,    -6e-7,
  997.         1.71548692e2, 1.4844328,    2.372e-4,    -6.1e-7,
  998.         4.63444e-2,    -2.658e-5,    7.7e-8,        0.0,
  999.         .772464,    6.253e-4,    3.95e-5,    0.0,
  1000.         73.477111,    .4986678,    1.3117e-3,    0.0,
  1001.         19.21814,    65.8,        -7.19
  1002.     },
  1003.  
  1004.     {   /*     neptune... */
  1005.  
  1006.         84.457994,    .6107942056,    3.205e-4,    -6e-7,
  1007.         4.6727364e1, 1.4245744,    3.9082e-4,    -6.05e-7,
  1008.         8.99704e-3,    6.33e-6,    -2e-9,        0.0,
  1009.         1.779242,    -9.5436e-3,    -9.1e-6,    0.0,
  1010.         130.681389,    1.098935,    2.4987e-4,    -4.718e-6,
  1011.         30.10957,    62.2,        -6.87
  1012.     },
  1013.  
  1014.     {   /*     pluto...(osculating 1984 jan 21) */
  1015.  
  1016.         95.3113544,    .3980332167,    0.0,        0.0,
  1017.         224.017,    0.0,        0.0,        0.0,
  1018.         .25515,    0.0,        0.0,        0.0,
  1019.         17.1329,    0.0,        0.0,        0.0,
  1020.         110.191,    0.0,        0.0,        0.0,
  1021.         39.8151,    8.2,        -1.0
  1022.     }
  1023. };
  1024.  
  1025. /* given a modified Julian date, mjd, return the elements for the mean orbit
  1026.  *   at that instant of all the major planets, together with their
  1027.  *   mean daily motions in longitude, angular diameter and standard visual
  1028.  *   magnitude.
  1029.  * plan[i][j] contains all the values for all the planets at mjd, such that
  1030.  *   i = 0..7: mercury, venus, mars, jupiter, saturn, unranus, neptune, pluto;
  1031.  *   j = 0..8: mean longitude, mean daily motion in longitude, longitude of 
  1032.  *     the perihelion, eccentricity, inclination, longitude of the ascending
  1033.  *     node, length of the semi-major axis, angular diameter from 1 AU, and
  1034.  *     the standard visual magnitude (see elements[][] comment, above).
  1035.  */
  1036. pelement (mjd, plan)
  1037. float mjd;
  1038. float plan[8][9];
  1039. {
  1040.     register float *ep, *pp;
  1041.     register float t = mjd/36525.;
  1042.     float aa;
  1043.     int planet, i;
  1044.  
  1045.     for (planet = 0; planet < 8; planet++) {
  1046.         ep = elements[planet];
  1047.         pp = plan[planet];
  1048.         aa = ep[1]*t;
  1049.         pp[0] = ep[0] + 360.*(aa-(int)aa) + (ep[3]*t + ep[2])*t*t;
  1050.         range (pp, 360.);
  1051.         pp[1] = (ep[1]*9.856263e-3) + (ep[2] + ep[3])/36525;
  1052.  
  1053.         for (i = 4; i < 20; i += 4)
  1054.         pp[i/4+1] = ((ep[i+3]*t + ep[i+2])*t + ep[i+1])*t + ep[i+0];
  1055.  
  1056.         pp[6] = ep[20];
  1057.         pp[7] = ep[21];
  1058.         pp[8] = ep[22];
  1059.     }
  1060. }
  1061. xXx
  1062. echo extracting plans.c
  1063. cat > plans.c << 'xXx'
  1064. #include <stdio.h>
  1065. #include <math.h>
  1066. #include "astro.h"
  1067.  
  1068. #define    TWOPI        (2*PI)
  1069. #define    mod2PI(x)    ((x) - (int)((x)/TWOPI)*TWOPI)
  1070.  
  1071. /* given a modified Julian date, mjd, and a planet, p, find:
  1072.  *   lpd0: heliocentric longitude, 
  1073.  *   psi0: heliocentric latitude,
  1074.  *   rp0:  distance from the sun to the planet, 
  1075.  *   rho0: distance from the Earth to the planet,
  1076.  *         none corrected for light time, ie, they are the true values for the
  1077.  *         given instant.
  1078.  *   lam:  geocentric ecliptic longitude, 
  1079.  *   bet:  geocentric ecliptic latitude,
  1080.  *         each corrected for light time, ie, they are the apparent values as
  1081.  *       seen from the center of the Earth for the given instant.
  1082.  *   dia:  angular diameter in arcsec at 1 AU, 
  1083.  *   mag:  visual magnitude when 1 AU from sun and earth at 0 phase angle.
  1084.  *
  1085.  * all angles are in radians, all distances in AU.
  1086.  * the mean orbital elements are found by calling pelement(), then mutual
  1087.  *   perturbation corrections are applied as necessary.
  1088.  *
  1089.  * corrections for nutation and abberation must be made by the caller. The RA 
  1090.  *   and DEC calculated from the fully-corrected ecliptic coordinates are then
  1091.  *   the apparent geocentric coordinates. Further corrections can be made, if
  1092.  *   required, for atmospheric refraction and geocentric parallax although the
  1093.  *   intrinsic error herein of about 10 arcseconds is usually the dominant
  1094.  *   error at this stage.
  1095.  * TODO: combine the several intermediate expressions when get a good compiler.
  1096.  */
  1097. plans (mjd, p, lpd0, psi0, rp0, rho0, lam, bet, dia, mag)
  1098. float mjd;
  1099. int p;
  1100. float *lpd0, *psi0, *rp0, *rho0, *lam, *bet, *dia, *mag;
  1101. {
  1102.     static float plan[8][9];
  1103.     static lastmjd = -10000;
  1104.     float dl;    /* perturbation correction for longitude */
  1105.     float dr;    /*  "   orbital radius */
  1106.     float dml;    /*  "   mean longitude */
  1107.     float ds;    /*  "   eccentricity */
  1108.     float dm;    /*  "   mean anomaly */
  1109.     float da;    /*  "   semi-major axis */
  1110.     float dhl;    /*  "   heliocentric longitude */
  1111.     float lsn, rsn;    /* true geocentric longitude of sun and sun-earth rad */
  1112.     float mas;    /* mean anomaly of the sun */
  1113.     float re;    /* radius of earth's orbit */
  1114.     float lg;    /* longitude of earth */
  1115.     float map[8];    /* array of mean anomalies for each planet */
  1116.     float lpd, psi, rp, rho;
  1117.     float ll, sll, cll;
  1118.     float t;
  1119.     float dt;
  1120.     int pass;
  1121.     int j;
  1122.     float s, ma;
  1123.     float nu, ea;
  1124.     float lp, om;
  1125.     float lo, slo, clo;
  1126.     float inc, y;
  1127.     float spsi, cpsi;
  1128.     float rpd;
  1129.  
  1130.     /* only need to fill in plan[] once for a given mjd */
  1131.     if (mjd != lastmjd) {
  1132.         pelement (mjd, plan);
  1133.         lastmjd = mjd;
  1134.     }
  1135.  
  1136.     dt = 0;
  1137.     t = mjd/36525.;
  1138.     sun (mjd, &lsn, &rsn);
  1139.     masun (mjd, &mas);
  1140.         re = rsn;
  1141.     lg = lsn+PI;
  1142.  
  1143.     /* first find the true position of the planet at mjd.
  1144.      * then repeat a second time for a slightly different time based
  1145.      * on the position found in the first pass to account for light-travel
  1146.      * time.
  1147.      */
  1148.     for (pass = 0; pass < 2; pass++) {
  1149.  
  1150.         for (j = 0; j < 8; j++)
  1151.         map[j] = degrad(plan[j][0]-plan[j][2]-dt*plan[j][1]);
  1152.  
  1153.         /* set initial corrections to 0.
  1154.          * then modify as necessary for the planet of interest.
  1155.          */
  1156.         dl = 0;
  1157.         dr = 0;
  1158.         dml = 0;
  1159.         ds = 0;
  1160.         dm = 0;
  1161.         da = 0;
  1162.         dhl = 0;
  1163.  
  1164.         switch (p) {
  1165.  
  1166.         case MERCURY:
  1167.         p_mercury (map, &dl, &dr);
  1168.         break;
  1169.  
  1170.         case VENUS:
  1171.         p_venus (t, mas, map, &dl, &dr, &dml, &dm);
  1172.         break;
  1173.  
  1174.         case MARS:
  1175.         p_mars (mas, map, &dl, &dr, &dml, &dm);
  1176.         break;
  1177.  
  1178.         case JUPITER:
  1179.         p_jupiter (t, plan[p][3], &dml, &ds, &dm, &da);
  1180.         break;
  1181.  
  1182.         case SATURN:
  1183.         p_saturn (t, plan[p][3], &dml, &ds, &dm, &da, &dhl);
  1184.         break;
  1185.  
  1186.         case URANUS:
  1187.         p_uranus (t, plan[p][3], &dl, &dr, &dml, &ds, &dm, &da, &dhl);
  1188.         break;
  1189.  
  1190.         case NEPTUNE:
  1191.         p_neptune (t, plan[p][3], &dl, &dr, &dml, &ds, &dm, &da, &dhl);
  1192.         break;
  1193.  
  1194.         case PLUTO:
  1195.         /* no perturbation theory for pluto */
  1196.         break;
  1197.         }
  1198.  
  1199.         s = plan[p][3]+ds;
  1200.         ma = map[p]+dm;
  1201.         anomaly (ma, s, &nu, &ea);
  1202.         rp = (plan[p][6]+da)*(1-s*s)/(1+s*cos(nu));
  1203.         lp = raddeg(nu)+plan[p][2]+raddeg(dml-dm);
  1204.         lp = degrad(lp);
  1205.         om = degrad(plan[p][5]);
  1206.         lo = lp-om;
  1207.         slo = sin(lo);
  1208.         clo = cos(lo);
  1209.         inc = degrad(plan[p][4]);
  1210.         rp = rp+dr;
  1211.         spsi = slo*sin(inc);
  1212.         y = slo*cos(inc);
  1213.         psi = asin(spsi)+dhl;
  1214.         spsi = sin(psi);
  1215.         lpd = atan(y/clo)+om+degrad(dl);
  1216.         if (clo<0) lpd += PI;
  1217.         if (lpd>TWOPI) lpd -= TWOPI;
  1218.         cpsi = cos(psi);
  1219.         rpd = rp*cpsi;
  1220.         ll = lpd-lg;
  1221.         rho = sqrt(re*re+rp*rp-2*re*rp*cpsi*cos(ll));
  1222.  
  1223.         /* when we view a planet we see it in the position it occupied
  1224.          * dt hours ago, where rho is the distance between it and earth,
  1225.          * in AU. use this as the new time for the next pass.
  1226.          */
  1227.         dt = rho*5.775518e-3;
  1228.  
  1229.         if (pass == 0) {
  1230.         /* save heliocentric coordinates after first pass since, being
  1231.          * true, they are NOT to be corrected for light-travel time.
  1232.          */
  1233.         *lpd0 = lpd;
  1234.         *psi0 = psi;
  1235.         *rp0 = rp;
  1236.         *rho0 = rho;
  1237.         }
  1238.     }
  1239.  
  1240.         sll = sin(ll);
  1241.     cll = cos(ll);
  1242.         if (p < MARS) 
  1243.         *lam = atan(-1*rpd*sll/(re-rpd*cll))+lg+PI;
  1244.     else
  1245.         *lam = atan(re*sll/(rpd-re*cll))+lpd;
  1246.     range (lam, TWOPI);
  1247.         *bet = atan(rpd*spsi*sin(*lam-lpd)/(cpsi*re*sll));
  1248.     *dia = plan[p][7];
  1249.     *mag = plan[p][8];
  1250. }
  1251.  
  1252. /* set auxilliary variables used for jupiter, saturn, uranus, and neptune */
  1253. static
  1254. aux_jsun (t, j1, j2, j3, j4, j5, j6)
  1255. float t;
  1256. float *j1, *j2, *j3, *j4, *j5, *j6;
  1257. {
  1258.         *j1 = t/5+0.1;
  1259.         *j2 = mod2PI(4.14473+5.29691e1*t);
  1260.         *j3 = mod2PI(4.641118+2.132991e1*t);
  1261.         *j4 = mod2PI(4.250177+7.478172*t);
  1262.         *j5 = 5 * *j3 - 2 * *j2;
  1263.     *j6 = 2 * *j2 - 6 * *j3 + 3 * *j4;
  1264. }
  1265.  
  1266. /* find the mean anomaly of the sun at mjd.
  1267.  * this is the same as that used in sun() but when it was converted to C it
  1268.  * was not known it would be required outside that routine.
  1269.  * TODO: add an argument to sun() to return mas and eliminate this routine.
  1270.  */
  1271. static
  1272. masun (mjd, mas)
  1273. float mjd;
  1274. float *mas;
  1275. {
  1276.     float t, t2;
  1277.     float a, b;
  1278.  
  1279.     t = mjd/36525;
  1280.     t2 = t*t;
  1281.     a = 9.999736042e1*t;
  1282.     b = 360.*(a-(int)a);
  1283.     *mas = degrad (3.5847583e2-(1.5e-4+3.3e-6*t)*t2+b);
  1284. }
  1285.  
  1286. /* perturbations for mercury */
  1287. static
  1288. p_mercury (map, dl, dr)
  1289. float map[];
  1290. float *dl, *dr;
  1291. {
  1292.     *dl = 2.04e-3*cos(5*map[2-1]-2*map[1-1]+2.1328e-1)+
  1293.          1.03e-3*cos(2*map[2-1]-map[1-1]-2.8046)+
  1294.          9.1e-4*cos(2*map[3]-map[1-1]-6.4582e-1)+
  1295.          7.8e-4*cos(5*map[2-1]-3*map[1-1]+1.7692e-1);
  1296.  
  1297.     *dr = 7.525e-6*cos(2*map[3]-map[1-1]+9.25251e-1)+
  1298.          6.802e-6*cos(5*map[2-1]-3*map[1-1]-4.53642)+
  1299.          5.457e-6*cos(2*map[2-1]-2*map[1-1]-1.24246)+
  1300.          3.569e-6*cos(5*map[2-1]-map[1-1]-1.35699);
  1301. }
  1302.  
  1303. /* ....venus */
  1304. static
  1305. p_venus (t, mas, map, dl, dr, dml, dm)
  1306. float t, mas, map[];
  1307. float *dl, *dr, *dml, *dm;
  1308. {
  1309.     *dml = degrad (7.7e-4*sin(4.1406+t*2.6227));
  1310.     *dm = *dml;
  1311.  
  1312.     *dl = 3.13e-3*cos(2*mas-2*map[2-1]-2.587)+
  1313.          1.98e-3*cos(3*mas-3*map[2-1]+4.4768e-2)+
  1314.          1.36e-3*cos(mas-map[2-1]-2.0788)+
  1315.          9.6e-4*cos(3*mas-2*map[2-1]-2.3721)+
  1316.          8.2e-4*cos(map[3]-map[2-1]-3.6318);
  1317.  
  1318.     *dr = 2.2501e-5*cos(2*mas-2*map[2-1]-1.01592)+
  1319.          1.9045e-5*cos(3*mas-3*map[2-1]+1.61577)+
  1320.          6.887e-6*cos(map[3]-map[2-1]-2.06106)+
  1321.          5.172e-6*cos(mas-map[2-1]-5.08065e-1)+
  1322.          3.62e-6*cos(5*mas-4*map[2-1]-1.81877)+
  1323.          3.283e-6*cos(4*mas-4*map[2-1]+1.10851)+
  1324.          3.074e-6*cos(2*map[3]-2*map[2-1]-9.62846e-1);
  1325. }
  1326.  
  1327. /* ....mars */
  1328. static
  1329. p_mars (mas, map, dl, dr, dml, dm)
  1330. float mas, map[];
  1331. float *dl, *dr, *dml, *dm;
  1332. {
  1333.     float a;
  1334.  
  1335.     a = 3*map[3]-8*map[2]+4*mas;
  1336.     *dml = degrad (-1*(1.133e-2*sin(a)+9.33e-3*cos(a)));
  1337.     *dm = *dml;
  1338.  
  1339.     *dl = 7.05e-3*cos(map[3]-map[2]-8.5448e-1)+
  1340.          6.07e-3*cos(2*map[3]-map[2]-3.2873)+
  1341.          4.45e-3*cos(2*map[3]-2*map[2]-3.3492)+
  1342.          3.88e-3*cos(mas-2*map[2]+3.5771e-1)+
  1343.          2.38e-3*cos(mas-map[2]+6.1256e-1)+
  1344.          2.04e-3*cos(2*mas-3*map[2]+2.7688)+
  1345.          1.77e-3*cos(3*map[2]-map[2-1]-1.0053)+
  1346.          1.36e-3*cos(2*mas-4*map[2]+2.6894)+
  1347.          1.04e-3*cos(map[3]+3.0749e-1);
  1348.  
  1349.     *dr = 5.3227e-5*cos(map[3]-map[2]+7.17864e-1)+
  1350.          5.0989e-5*cos(2*map[3]-2*map[2]-1.77997)+
  1351.          3.8278e-5*cos(2*map[3]-map[2]-1.71617)+
  1352.          1.5996e-5*cos(mas-map[2]-9.69618e-1)+
  1353.          1.4764e-5*cos(2*mas-3*map[2]+1.19768)+
  1354.          8.966e-6*cos(map[3]-2*map[2]+7.61225e-1)+
  1355.          7.914e-6*cos(3*map[3]-2*map[2]-2.43887)+
  1356.          7.004e-6*cos(2*map[3]-3*map[2]-1.79573)+
  1357.          6.62e-6*cos(mas-2*map[2]+1.97575)+
  1358.          4.93e-6*cos(3*map[3]-3*map[2]-1.33069)+
  1359.          4.693e-6*cos(3*mas-5*map[2]+3.32665)+
  1360.          4.571e-6*cos(2*mas-4*map[2]+4.27086)+
  1361.          4.409e-6*cos(3*map[3]-map[2]-2.02158);
  1362. }
  1363.  
  1364. /* ....jupiter */
  1365. static
  1366. p_jupiter (t, s, dml, ds, dm, da)
  1367. float t, s;
  1368. float *dml, *ds, *dm, *da;
  1369. {
  1370.     float dp;
  1371.     float j1, j2, j3, j4, j5, j6, j7;
  1372.     float sj3, cj3, s2j3, c2j3;
  1373.         float sj5, cj5, s2j5;
  1374.     float sj6;
  1375.         float sj7, cj7, s2j7, c2j7, s3j7, c3j7, s4j7, c4j7, c5j7;
  1376.  
  1377.     aux_jsun (t, &j1, &j2, &j3, &j4, &j5, &j6);
  1378.         j7 = j3-j2;
  1379.     sj3 = sin(j3);
  1380.     cj3 = cos(j3);
  1381.         s2j3 = sin(2*j3);
  1382.     c2j3 = cos(2*j3);
  1383.         sj5 = sin(j5);
  1384.     cj5 = cos(j5);
  1385.         s2j5 = sin(2*j5);
  1386.     sj6 = sin(j6);
  1387.         sj7 = sin(j7);
  1388.     cj7 = cos(j7);
  1389.         s2j7 = sin(2*j7);
  1390.     c2j7 = cos(2*j7);
  1391.         s3j7 = sin(3*j7);
  1392.     c3j7 = cos(3*j7);
  1393.         s4j7 = sin(4*j7);
  1394.     c4j7 = cos(4*j7);
  1395.         c5j7 = cos(5*j7);
  1396.  
  1397.     *dml = (3.31364e-1-(1.0281e-2+4.692e-3*j1)*j1)*sj5+
  1398.           (3.228e-3-(6.4436e-2-2.075e-3*j1)*j1)*cj5-
  1399.           (3.083e-3+(2.75e-4-4.89e-4*j1)*j1)*s2j5+
  1400.           2.472e-3*sj6+1.3619e-2*sj7+1.8472e-2*s2j7+6.717e-3*s3j7+
  1401.           2.775e-3*s4j7+6.417e-3*s2j7*sj3+
  1402.           (7.275e-3-1.253e-3*j1)*sj7*sj3+
  1403.           2.439e-3*s3j7*sj3-(3.5681e-2+1.208e-3*j1)*sj7*cj3-
  1404.           3.767e-3*c2j7*sj3-(3.3839e-2+1.125e-3*j1)*cj7*sj3-
  1405.           4.261e-3*s2j7*cj3+
  1406.           (1.161e-3*j1-6.333e-3)*cj7*cj3+
  1407.           2.178e-3*cj3-6.675e-3*c2j7*cj3-2.664e-3*c3j7*cj3-
  1408.           2.572e-3*sj7*s2j3-3.567e-3*s2j7*s2j3+2.094e-3*cj7*c2j3+
  1409.           3.342e-3*c2j7*c2j3;
  1410.     *dml = degrad(*dml);
  1411.  
  1412.     *ds = (3606+(130-43*j1)*j1)*sj5+(1289-580*j1)*cj5-6764*sj7*sj3-
  1413.          1110*s2j7*sj3-224*s3j7*sj3-204*sj3+(1284+116*j1)*cj7*sj3+
  1414.          188*c2j7*sj3+(1460+130*j1)*sj7*cj3+224*s2j7*cj3-817*cj3+
  1415.          6074*cj3*cj7+992*c2j7*cj3+
  1416.          508*c3j7*cj3+230*c4j7*cj3+108*c5j7*cj3-
  1417.          (956+73*j1)*sj7*s2j3+448*s2j7*s2j3+137*s3j7*s2j3+
  1418.          (108*j1-997)*cj7*s2j3+480*c2j7*s2j3+148*c3j7*s2j3+
  1419.          (99*j1-956)*sj7*c2j3+490*s2j7*c2j3+
  1420.          158*s3j7*c2j3+179*c2j3+(1024+75*j1)*cj7*c2j3-
  1421.          437*c2j7*c2j3-132*c3j7*c2j3;
  1422.     *ds *= 1e-7;
  1423.  
  1424.     dp = (7.192e-3-3.147e-3*j1)*sj5-4.344e-3*sj3+
  1425.          (j1*(1.97e-4*j1-6.75e-4)-2.0428e-2)*cj5+
  1426.          3.4036e-2*cj7*sj3+(7.269e-3+6.72e-4*j1)*sj7*sj3+
  1427.          5.614e-3*c2j7*sj3+2.964e-3*c3j7*sj3+3.7761e-2*sj7*cj3+
  1428.          6.158e-3*s2j7*cj3-
  1429.          6.603e-3*cj7*cj3-5.356e-3*sj7*s2j3+2.722e-3*s2j7*s2j3+
  1430.          4.483e-3*cj7*s2j3-2.642e-3*c2j7*s2j3+4.403e-3*sj7*c2j3-
  1431.          2.536e-3*s2j7*c2j3+5.547e-3*cj7*c2j3-2.689e-3*c2j7*c2j3;
  1432.  
  1433.     *dm = *dml-(degrad(dp)/s);
  1434.  
  1435.     *da = 205*cj7-263*cj5+693*c2j7+312*c3j7+147*c4j7+299*sj7*sj3+
  1436.          181*c2j7*sj3+204*s2j7*cj3+111*s3j7*cj3-337*cj7*cj3-
  1437.          111*c2j7*cj3;
  1438.     *da *= 1e-6;
  1439. }
  1440.  
  1441. /* ....saturn */
  1442. static
  1443. p_saturn (t, s, dml, ds, dm, da, dhl)
  1444. float t, s;
  1445. float *dml, *ds, *dm, *da, *dhl;
  1446. {
  1447.     float dp;
  1448.     float j1, j2, j3, j4, j5, j6, j7, j8;
  1449.     float sj3, cj3, s2j3, c2j3, s3j3, c3j3, s4j3, c4j3;
  1450.         float sj5, cj5, s2j5, c2j5;
  1451.     float sj6;
  1452.         float sj7, cj7, s2j7, c2j7, s3j7, c3j7, s4j7, c4j7, c5j7, s5j7;
  1453.     float s2j8, c2j8, s3j8, c3j8;
  1454.  
  1455.     aux_jsun (t, &j1, &j2, &j3, &j4, &j5, &j6);
  1456.         j7 = j3-j2;
  1457.     sj3 = sin(j3);
  1458.     cj3 = cos(j3);
  1459.         s2j3 = sin(2*j3);
  1460.     c2j3 = cos(2*j3);
  1461.         sj5 = sin(j5);
  1462.     cj5 = cos(j5);
  1463.         s2j5 = sin(2*j5);
  1464.     sj6 = sin(j6);
  1465.         sj7 = sin(j7);
  1466.     cj7 = cos(j7);
  1467.         s2j7 = sin(2*j7);
  1468.     c2j7 = cos(2*j7);
  1469.         s3j7 = sin(3*j7);
  1470.     c3j7 = cos(3*j7);
  1471.         s4j7 = sin(4*j7);
  1472.     c4j7 = cos(4*j7);
  1473.         c5j7 = cos(5*j7);
  1474.  
  1475.     s3j3 = sin(3*j3);
  1476.     c3j3 = cos(3*j3);
  1477.     s4j3 = sin(4*j3);
  1478.     c4j3 = cos(4*j3);
  1479.     c2j5 = cos(2*j5);
  1480.     s5j7 = sin(5*j7);
  1481.     j8 = j4-j3;
  1482.     s2j8 = sin(2*j8);
  1483.     c2j8 = cos(2*j8);
  1484.     s3j8 = sin(3*j8);
  1485.     c3j8 = cos(3*j8);
  1486.  
  1487.     *dml = 7.581e-3*s2j5-7.986e-3*sj6-1.48811e-1*sj7-4.0786e-2*s2j7-
  1488.           (8.14181e-1-(1.815e-2-1.6714e-2*j1)*j1)*sj5-
  1489.           (1.0497e-2-(1.60906e-1-4.1e-3*j1)*j1)*cj5-1.5208e-2*s3j7-
  1490.           6.339e-3*s4j7-6.244e-3*sj3-1.65e-2*s2j7*sj3+
  1491.           (8.931e-3+2.728e-3*j1)*sj7*sj3-5.775e-3*s3j7*sj3+
  1492.           (8.1344e-2+3.206e-3*j1)*cj7*sj3+1.5019e-2*c2j7*sj3+
  1493.           (8.5581e-2+2.494e-3*j1)*sj7*cj3+1.4394e-2*c2j7*cj3+
  1494.           (2.5328e-2-3.117e-3*j1)*cj7*cj3+
  1495.           6.319e-3*c3j7*cj3+6.369e-3*sj7*s2j3+9.156e-3*s2j7*s2j3+
  1496.           7.525e-3*s3j8*s2j3-5.236e-3*cj7*c2j3-7.736e-3*c2j7*c2j3-
  1497.           7.528e-3*c3j8*c2j3;
  1498.     *dml = degrad(*dml);
  1499.  
  1500.     *ds = (-7927+(2548+91*j1)*j1)*sj5+(13381+(1226-253*j1)*j1)*cj5+
  1501.          (248-121*j1)*s2j5-(305+91*j1)*c2j5+412*s2j7+12415*sj3+
  1502.          (390-617*j1)*sj7*sj3+(165-204*j1)*s2j7*sj3+26599*cj7*sj3-
  1503.          4687*c2j7*sj3-1870*c3j7*sj3-821*c4j7*sj3-
  1504.          377*c5j7*sj3+497*c2j8*sj3+(163-611*j1)*cj3-
  1505.          12696*sj7*cj3-4200*s2j7*cj3-1503*s3j7*cj3-619*s4j7*cj3-
  1506.          268*s5j7*cj3-(282+1306*j1)*cj7*cj3+(-86+230*j1)*c2j7*cj3+
  1507.          461*s2j8*cj3-350*s2j3+(2211-286*j1)*sj7*s2j3-
  1508.          2208*s2j7*s2j3-568*s3j7*s2j3-346*s4j7*s2j3-
  1509.          (2780+222*j1)*cj7*s2j3+(2022+263*j1)*c2j7*s2j3+248*c3j7*s2j3+
  1510.          242*s3j8*s2j3+467*c3j8*s2j3-490*c2j3-(2842+279*j1)*sj7*c2j3+
  1511.          (128+226*j1)*s2j7*c2j3+224*s3j7*c2j3+
  1512.          (-1594+282*j1)*cj7*c2j3+(2162-207*j1)*c2j7*c2j3+
  1513.          561*c3j7*c2j3+343*c4j7*c2j3+469*s3j8*c2j3-242*c3j8*c2j3-
  1514.          205*sj7*s3j3+262*s3j7*s3j3+208*cj7*c3j3-271*c3j7*c3j3-
  1515.          382*c3j7*s4j3-376*s3j7*c4j3;
  1516.     *ds *= 1e-7;
  1517.  
  1518.     dp = (7.7108e-2+(7.186e-3-1.533e-3*j1)*j1)*sj5-7.075e-3*sj7+
  1519.          (4.5803e-2-(1.4766e-2+5.36e-4*j1)*j1)*cj5-7.2586e-2*cj3-
  1520.          7.5825e-2*sj7*sj3-2.4839e-2*s2j7*sj3-8.631e-3*s3j7*sj3-
  1521.          1.50383e-1*cj7*cj3+2.6897e-2*c2j7*cj3+1.0053e-2*c3j7*cj3-
  1522.          (1.3597e-2+1.719e-3*j1)*sj7*s2j3+1.1981e-2*s2j7*c2j3-
  1523.          (7.742e-3-1.517e-3*j1)*cj7*s2j3+
  1524.          (1.3586e-2-1.375e-3*j1)*c2j7*c2j3-
  1525.          (1.3667e-2-1.239e-3*j1)*sj7*c2j3+
  1526.          (1.4861e-2+1.136e-3*j1)*cj7*c2j3-
  1527.          (1.3064e-2+1.628e-3*j1)*c2j7*c2j3;
  1528.  
  1529.     *dm = *dml-(degrad(dp)/s);
  1530.  
  1531.     *da = 572*sj5-1590*s2j7*cj3+2933*cj5-647*s3j7*cj3+33629*cj7-
  1532.          344*s4j7*cj3-3081*c2j7+2885*cj7*cj3-1423*c3j7+
  1533.          (2172+102*j1)*c2j7*cj3-671*c4j7+296*c3j7*cj3-320*c5j7-
  1534.          267*s2j7*s2j3+1098*sj3-778*cj7*s2j3-2812*sj7*sj3+
  1535.          495*c2j7*s2j3+688*s2j7*sj3+250*c3j7*s2j3-393*s3j7*sj3-
  1536.          856*sj7*c2j3-228*s4j7*sj3+441*s2j7*c2j3+2138*cj7*sj3+
  1537.          296*c2j7*c2j3-999*c2j7*sj3+211*c3j7*c2j3-642*c3j7*sj3-
  1538.          427*sj7*s3j3-325*c4j7*sj3+398*s3j7*s3j3-890*cj3+
  1539.          344*cj7*c3j3+2206*sj7*cj3-427*c3j7*c3j3;
  1540.     *da *= 1e-6;
  1541.  
  1542.     *dhl = 7.47e-4*cj7*sj3+1.069e-3*cj7*cj3+2.108e-3*s2j7*s2j3+
  1543.           1.261e-3*c2j7*s2j3+1.236e-3*s2j7*c2j3-2.075e-3*c2j7*c2j3;
  1544.     *dhl = degrad(*dhl);
  1545. }
  1546.  
  1547. /* ....uranus */
  1548. static
  1549. p_uranus (t, s, dl, dr, dml, ds, dm, da, dhl)
  1550. float t, s;
  1551. float *dl, *dr, *dml, *ds, *dm, *da, *dhl;
  1552. {
  1553.     float dp;
  1554.     float j1, j2, j3, j4, j5, j6;
  1555.     float j8, j9, j10, j11, j12;
  1556.     float sj4, cj4, s2j4, c2j4;
  1557.     float sj9, cj9, s2j9, c2j9;
  1558.     float sj11, cj11;
  1559.  
  1560.     aux_jsun (t, &j1, &j2, &j3, &j4, &j5, &j6);
  1561.  
  1562.         j8 = mod2PI(1.46205+3.81337*t);
  1563.         j9 = 2*j8-j4;
  1564.     sj9 = sin(j9);
  1565.     cj9 = cos(j9);
  1566.         s2j9 = sin(2*j9);
  1567.     c2j9 = cos(2*j9);
  1568.  
  1569.     j10 = j4-j2;
  1570.     j11 = j4-j3;
  1571.     j12 = j8-j4;
  1572.  
  1573.     *dml = (8.64319e-1-1.583e-3*j1)*sj9+(8.2222e-2-6.833e-3*j1)*cj9+
  1574.           3.6017e-2*s2j9-3.019e-3*c2j9+8.122e-3*sin(j6);
  1575.     *dml = degrad(*dml);
  1576.  
  1577.     dp = 1.20303e-1*sj9+6.197e-3*s2j9+(1.9472e-2-9.47e-4*j1)*cj9;
  1578.     *dm = *dml-(degrad(dp)/s);
  1579.  
  1580.     *ds = (163*j1-3349)*sj9+20981*cj9+1311*c2j9;
  1581.     *ds *= 1e-7;
  1582.  
  1583.     *da = -3.825e-3*cj9;
  1584.  
  1585.     *dl = (1.0122e-2-9.88e-4*j1)*sin(j4+j11)+
  1586.          (-3.8581e-2+(2.031e-3-1.91e-3*j1)*j1)*cos(j4+j11)+
  1587.          (3.4964e-2-(1.038e-3-8.68e-4*j1)*j1)*cos(2*j4+j11)+
  1588.          5.594e-3*sin(j4+3*j12)-1.4808e-2*sin(j10)-
  1589.          5.794e-3*sin(j11)+2.347e-3*cos(j11)+9.872e-3*sin(j12)+
  1590.          8.803e-3*sin(2*j12)-4.308e-3*sin(3*j12);
  1591.  
  1592.     sj11 = sin(j11);
  1593.     cj11 = cos(j11);
  1594.     sj4 = sin(j4);
  1595.     cj4 = cos(j4);
  1596.     s2j4 = sin(2*j4);
  1597.     c2j4 = cos(2*j4);
  1598.     *dhl = (4.58e-4*sj11-6.42e-4*cj11-5.17e-4*cos(4*j12))*sj4-
  1599.           (3.47e-4*sj11+8.53e-4*cj11+5.17e-4*sin(4*j11))*cj4+
  1600.           4.03e-4*(cos(2*j12)*s2j4+sin(2*j12)*c2j4);
  1601.     *dhl = degrad(*dhl);
  1602.  
  1603.     *dr = -25948+4985*cos(j10)-1230*cj4+3354*cos(j11)+904*cos(2*j12)+
  1604.          894*(cos(j12)-cos(3*j12))+(5795*cj4-1165*sj4+1388*c2j4)*sj11+
  1605.          (1351*cj4+5702*sj4+1388*s2j4)*cos(j11);
  1606.     *dr *= 1e-6;
  1607. }
  1608.  
  1609. /* ....neptune */
  1610. static
  1611. p_neptune (t, s, dl, dr, dml, ds, dm, da, dhl)
  1612. float t, s;
  1613. float *dl, *dr, *dml, *ds, *dm, *da, *dhl;
  1614. {
  1615.     float dp;
  1616.     float j1, j2, j3, j4, j5, j6;
  1617.     float j8, j9, j10, j11, j12;
  1618.     float sj8, cj8;
  1619.     float sj9, cj9, s2j9, c2j9;
  1620.     float s2j12, c2j12;
  1621.  
  1622.     aux_jsun (t, &j1, &j2, &j3, &j4, &j5, &j6);
  1623.  
  1624.         j8 = mod2PI(1.46205+3.81337*t);
  1625.         j9 = 2*j8-j4;
  1626.     sj9 = sin(j9);
  1627.     cj9 = cos(j9);
  1628.         s2j9 = sin(2*j9);
  1629.     c2j9 = cos(2*j9);
  1630.  
  1631.     j10 = j8-j2;
  1632.     j11 = j8-j3;
  1633.     j12 = j8-j4;
  1634.  
  1635.     *dml = (1.089e-3*j1-5.89833e-1)*sj9+(4.658e-3*j1-5.6094e-2)*cj9-
  1636.           2.4286e-2*s2j9;
  1637.     *dml = degrad(*dml);
  1638.  
  1639.     dp = 2.4039e-2*sj9-2.5303e-2*cj9+6.206e-3*s2j9-5.992e-3*c2j9;
  1640.  
  1641.     *dm = *dml-(degrad(dp)/s);
  1642.  
  1643.     *ds = 4389*sj9+1129*s2j9+4262*cj9+1089*c2j9;
  1644.     *ds *= 1e-7;
  1645.  
  1646.     *da = 8189*cj9-817*sj9+781*c2j9;
  1647.     *da *= 1e-6;
  1648.  
  1649.     s2j12 = sin(2*j12);
  1650.     c2j12 = cos(2*j12);
  1651.     sj8 = sin(j8);
  1652.     cj8 = cos(j8);
  1653.     *dl = -9.556e-3*sin(j10)-5.178e-3*sin(j11)+2.572e-3*s2j12-
  1654.          2.972e-3*c2j12*sj8-2.833e-3*s2j12*cj8;
  1655.  
  1656.     *dhl = 3.36e-4*c2j12*sj8+3.64e-4*s2j12*cj8;
  1657.     *dhl = degrad(*dhl);
  1658.  
  1659.     *dr = -40596+4992*cos(j10)+2744*cos(j11)+2044*cos(j12)+1051*c2j12;
  1660.     *dr *= 1e-6;
  1661. }
  1662. xXx
  1663. echo extracting precess.c
  1664. cat > precess.c << 'xXx'
  1665. #include <stdio.h>
  1666. #include <math.h>
  1667. #include "astro.h"
  1668.  
  1669. /* corrects ra and dec, both in radians, for precession from epoch 1 to epoch 2.
  1670.  * the epochs are given by their modified JDs, mjd1 and mjd2, respectively.
  1671.  * N.B. ra and dec are modifed IN PLACE.
  1672.  */
  1673. precess (mjd1, mjd2, ra, dec)
  1674. float mjd1, mjd2;    /* initial and final epoch modified JDs */
  1675. float *ra, *dec;    /* ra/dec for mjd1 in, for mjd2 out */
  1676. {
  1677.     static float lastmjd1, lastmjd2;
  1678.     static float m, n, nyrs;
  1679.     float dra, ddec;    /* ra and dec corrections */
  1680.  
  1681.     if (mjd1 != lastmjd1 || mjd2 != lastmjd2) {
  1682.         float t1, t2; /* Julian centuries of 36525 days since Jan .5 1900 */
  1683.         float m1, n1; /* "constants" for t1 */
  1684.         float m2, n2; /* "constants" for t2 */
  1685.         t1 = mjd1/36525.;
  1686.         t2 = mjd2/36525.;
  1687.         m1 = 3.07234+(1.86e-3*t1);
  1688.         n1 = 20.0468-(8.5e-3*t1);
  1689.         m2 = 3.07234+(1.86e-3*t2);
  1690.         n2 = 20.0468-(8.5e-3*t2);
  1691.         m = (m1+m2)/2;    /* average m for the two epochs */
  1692.         n = (n1+n2)/2;    /* average n for the two epochs */
  1693.         nyrs = (mjd2-mjd1)/365.2425;
  1694.         lastmjd1 = mjd1;
  1695.         lastmjd2 = mjd2;
  1696.     }
  1697.  
  1698.     dra = (m+(n*sin(*ra)*tan(*dec)/15))*7.272205e-5*nyrs;
  1699.     ddec = n*cos(*ra)*4.848137e-6*nyrs;
  1700.     *ra += dra;
  1701.     *dec += ddec;
  1702.     range (ra, 2*PI);
  1703. }
  1704. xXx
  1705. echo extracting refract.c
  1706. cat > refract.c << 'xXx'
  1707. #include <stdio.h>
  1708. #include <math.h>
  1709. #include "astro.h"
  1710.  
  1711. /* correct the true altitude, ta, for refraction to the apparent altitude, aa,
  1712.  * each in radians, given the local atmospheric pressure, pr, in mbars, and
  1713.  * the temperature, tr, in degrees C.
  1714.  */
  1715. refract (pr, tr, ta, aa)
  1716. float pr, tr;
  1717. float ta;
  1718. float *aa;
  1719. {
  1720.     float r;    /* refraction correction*/
  1721.  
  1722.         if (ta >= degrad(15.)) {
  1723.         /* model for altitudes at least 15 degrees above horizon */
  1724.             r = 7.888888e-5*pr/((273+tr)*tan(ta));
  1725.     } else if (ta > degrad(-5.)) {
  1726.         /* hairier model for altitudes at least -5 and below 15 degrees */
  1727.         float a, b, tadeg = raddeg(ta);
  1728.         a = ((2e-5*tadeg+1.96e-2)*tadeg+1.594e-1)*pr;
  1729.         b = (273+tr)*((8.45e-2*tadeg+5.05e-1)*tadeg+1);
  1730.         r = degrad(a/b);
  1731.     } else {
  1732.         /* do nothing if more than 5 degrees below horizon.
  1733.          */
  1734.         r = 0;
  1735.     }
  1736.  
  1737.     *aa  =  ta + r;
  1738. }
  1739.  
  1740. /* correct the apparent altitude, aa, for refraction to the true altitude, ta,
  1741.  * each in radians, given the local atmospheric pressure, pr, in mbars, and
  1742.  * the temperature, tr, in degrees C.
  1743.  */
  1744. unrefract (pr, tr, aa, ta)
  1745. float pr, tr;
  1746. float aa;
  1747. float *ta;
  1748. {
  1749.     float err;
  1750.     float appar;
  1751.     float true;
  1752.  
  1753.     /* iterative solution: search for the true that refracts to the
  1754.      *   given apparent.
  1755.      * since refract() is discontinuous at -5 degrees, there is a range
  1756.      *   of apparent altitudes between about -4.5 and -5 degrees that are
  1757.      *   not invertable (the graph of ap vs. true has a vertical step at
  1758.      *   true = -5). thus, the iteration just oscillates if it gets into
  1759.      *   this region. if this happens the iteration is forced to abort.
  1760.      *   of course, this makes unrefract() discontinuous too.
  1761.      */
  1762.     true = aa;
  1763.     do {
  1764.         refract (pr, tr, true, &appar);
  1765.         err = appar - aa;
  1766.         true -= err;
  1767.     } while (fabs(err) >= 1e-6 && true > degrad(-5));
  1768.  
  1769.     *ta = true;
  1770. }
  1771. xXx
  1772. echo extracting riset.c
  1773. cat > riset.c << 'xXx'
  1774. #include <stdio.h>
  1775. #include <math.h>
  1776. #include "astro.h"
  1777.  
  1778. /* given the true geocentric ra and dec of an object, in radians, the observer's
  1779.  *   latitude in radians, and a horizon displacement correction, dis, in radians
  1780.  *   find the local sidereal times and azimuths of rising and setting, lstr/s
  1781.  *   and azr/s, in radians, respectively.
  1782.  * dis is the vertical displacement from the true position of the horizon. it
  1783.  *   is positive if the apparent position is higher than the true position.
  1784.  *   said another way, it is positive if the shift causes the object to spend
  1785.  *   longer above the horizon. for example, atmospheric refraction is typically
  1786.  *   assumed to produce a vertical shift of 34 arc minutes at the horizon; dis
  1787.  *   would then take on the value +9.89e-3 (radians). On the other hand, if
  1788.  *   your horizon has hills such that your apparent horizon is, say, 1 degree
  1789.  *   above sea level, you would allow for this by setting dis to -1.75e-2
  1790.  *   (radians).
  1791.  * status: 0: normal; 1: never rises; -1: circumpolar; 2: trouble.
  1792.  */
  1793. riset (ra, dec, lat, dis, lstr, lsts, azr, azs, status)
  1794. float ra, dec;
  1795. float lat, dis;
  1796. float *lstr, *lsts;
  1797. float *azr, *azs;
  1798. int *status;
  1799. {
  1800.     static float lastlat = 0, slat = 0, clat = 1.0;
  1801.     float dec1, sdec, cdec, tdec;
  1802.     float psi, spsi, cpsi;
  1803.     float h, dh, ch;    /* hr angle, delta and cos */
  1804.  
  1805.     /* avoid needless sin/cos since latitude doesn't change often */
  1806.         if (lat != lastlat) {
  1807.         clat = cos(lat);
  1808.         slat = sin(lat);
  1809.         lastlat = lat;
  1810.     }
  1811.  
  1812.     /* can't cope with objects very near the celestial pole nor if we 
  1813.      * are located very near the earth's poles.
  1814.      */
  1815.     cdec = cos(dec);
  1816.         if (fabs(cdec*clat) < 1e-10) {
  1817.         /* trouble */
  1818.         *status = 2;
  1819.         return;
  1820.     }
  1821.  
  1822.         cpsi = slat/cdec;
  1823.         if (cpsi>1.0) cpsi = 1.0;
  1824.         else if (cpsi<-1.0) cpsi = -1.0;
  1825.         psi = acos(cpsi);
  1826.     spsi = sin(psi);
  1827.  
  1828.         dh = dis*spsi;
  1829.     dec1 = dec + (dis*cpsi);
  1830.         sdec = sin(dec1);
  1831.     tdec = tan(dec1);
  1832.  
  1833.         ch = (-1*slat*tdec)/clat;
  1834.  
  1835.         if (ch < -1) {
  1836.         /* circumpolar */
  1837.         *status = -1;
  1838.         return;
  1839.     }
  1840.         if (ch > 1) {
  1841.         /* never rises */
  1842.         *status = 1;
  1843.         return;
  1844.     }
  1845.  
  1846.         *status = 0;
  1847.     h = acos(ch)+dh;
  1848.  
  1849.         *lstr = 24+radhr(ra-h);
  1850.     *lsts = radhr(ra+h);
  1851.     range (lstr, 24.0);
  1852.     range (lsts, 24.0);
  1853.  
  1854.     *azr = acos(sdec/clat);
  1855.     range (azr, 2*PI);
  1856.         *azs = 2*PI - *azr;
  1857. }
  1858. xXx
  1859. echo extracting sex_dec.c
  1860. cat > sex_dec.c << 'xXx'
  1861. /* given hours (or degrees), hd, minutes, m, and seconds, s, 
  1862.  * return decimal hours (or degrees), *d.
  1863.  * in the case of hours (angles) < 0, only the first non-zero element should
  1864.  *   be negative.
  1865.  */
  1866. sex_dec (hd, m, s, d)
  1867. int hd, m, s;
  1868. float *d;
  1869. {
  1870.     int sign = 1;
  1871.  
  1872.     if (hd < 0) {
  1873.         sign = -1;
  1874.         hd = -hd;
  1875.     } else if (m < 0) {
  1876.         sign = -1;
  1877.         m = -m;
  1878.     } else if (s < 0) {
  1879.         sign = -1;
  1880.         s = -s;
  1881.     }
  1882.  
  1883.     *d = ((s/60.0 + m)/60.0 + hd) * sign;
  1884. }
  1885.  
  1886. /* given decimal hours (or degrees), d.
  1887.  * return nearest hours (or degrees), *hd, minutes, *m, and seconds, *s, 
  1888.  * each always non-negative; *isneg is set to 1 if d is < 0, else to 0.
  1889.  */
  1890. dec_sex (d, hd, m, s, isneg)
  1891. float d;
  1892. int *hd, *m, *s, *isneg;
  1893. {
  1894.     float min;
  1895.  
  1896.     if (d < 0) {
  1897.         *isneg = 1;
  1898.         d = -d;
  1899.     } else
  1900.         *isneg = 0;
  1901.  
  1902.     *hd = (int)d;
  1903.     min = (d - *hd)*60.;
  1904.     *m = (int)min;
  1905.     *s = (int)((min - *m)*60. + 0.5);
  1906.  
  1907.     if (*s == 60) {
  1908.         if ((*m += 1) == 60) {
  1909.         *hd += 1;
  1910.         *m = 0;
  1911.         }
  1912.         *s = 0;
  1913.     }
  1914.     /* no  negative 0's */
  1915.     if (*hd == 0 && *m == 0 && *s == 0)
  1916.         *isneg = 0;
  1917. }
  1918.  
  1919. /* insure 0 <= *v < r.
  1920.  */
  1921. range (v, r)
  1922. float *v, r;
  1923. {
  1924.     while (*v <  0) *v += r;
  1925.     while (*v >= r) *v -= r;
  1926. }
  1927. xXx
  1928. echo extracting sun.c
  1929. cat > sun.c << 'xXx'
  1930. #include <stdio.h>
  1931. #include <math.h>
  1932. #include "astro.h"
  1933.  
  1934. /* given the modified JD, mjd, return the true geocentric ecliptic longitude
  1935.  *   of the sun for the mean equinox of the date, *lsn, in radians, and the
  1936.  *   sun-earth distance, *rsn, in AU. (the true ecliptic latitude is never more
  1937.  *   than 1.2 arc seconds and so may be taken to be a constant 0.)
  1938.  * if the APPARENT ecliptic longitude is required, correct the longitude for
  1939.  *   nutation to the true equinox of date and for aberration (light travel time,
  1940.  *   approximately  -9.27e7/186000/(3600*24*365)*2*pi = -9.93e-5 radians).
  1941.  */
  1942. sun (mjd, lsn, rsn)
  1943. float mjd;
  1944. float *lsn, *rsn;
  1945. {
  1946.     float t, t2;
  1947.     float ls, ms;    /* mean longitude and mean anomoay */
  1948.     float s, nu, ea; /* eccentricity, true anomaly, eccentric anomaly */
  1949.     double a, b, a1, b1, c1, d1, e1, h1, dl, dr;
  1950.  
  1951.     t = mjd/36525.;
  1952.     t2 = t*t;
  1953.     a = 100.0021359*t;
  1954.     b = 360.*(a-(int)a);
  1955.     ls = 279.69668+.0003025*t2+b;
  1956.     a = 99.99736042000039*t;
  1957.     b = 360*(a-(int)a);
  1958.     ms = 358.47583-(.00015+.0000033*t)*t2+b;
  1959.     s = .016751-.0000418*t-1.26e-07*t2;
  1960.     anomaly (degrad(ms), s, &nu, &ea);
  1961.     a = 62.55209472000015*t;
  1962.     b = 360*(a-(int)a);
  1963.     a1 = degrad(153.23+b);
  1964.     a = 125.1041894*t;
  1965.     b = 360*(a-(int)a);
  1966.     b1 = degrad(216.57+b);
  1967.     a = 91.56766028*t;
  1968.     b = 360*(a = (int)a);
  1969.     c1 = degrad(312.69+b);
  1970.     a = 1236.853095*t;
  1971.     b = 360*(a-(int)a);
  1972.     d1 = degrad(350.74-.00144*t2+b);
  1973.     e1 = degrad(231.19+20.2*t);
  1974.     a = 183.1353208*t;
  1975.     b = 360*(a-(int)a);
  1976.     h1 = degrad(353.4+b);
  1977.     dl = .00134*cos(a1)+.00154*cos(b1)+.002*cos(c1)+.00179*sin(d1)+
  1978.                                 .00178*sin(e1);
  1979.     dr = 5.43e-06*sin(a1)+1.575e-05*sin(b1)+1.627e-05*sin(c1)+
  1980.                         3.076e-05*cos(d1)+9.27e-06*sin(h1);
  1981.     *lsn = nu+degrad(ls-ms+dl);
  1982.     *rsn = 1.0000002*(1-s*cos(ea))+dr;
  1983.     range (lsn, 2*PI);
  1984. }
  1985. xXx
  1986. echo extracting sunrs.c
  1987. cat > sunrs.c << 'xXx'
  1988. #include <stdio.h>
  1989. #include <math.h>
  1990. #include "astro.h"
  1991.  
  1992. /* given the modified julian date, mjd, geographical latitude, lat, and
  1993.  *   longitude, lng, each in radians with west lng < 0, and a horizon
  1994.  *   displacement, dis, find the utc times of sunrise, utcr, and sunset, utcs,
  1995.  *   and the azimuths of sunrise, azr, and sunset, azs. all angles in radians;
  1996.  *   azimuths are E of N.
  1997.  * times are for the upper limb, including effects of nutation and aberration.
  1998.  * displacement is additional amount added to suns altitude compared to its
  1999.  *   true location; see riset.c for more discussion. it can be used to
  2000.  *   account for irregular horizons, various refraction models, even times
  2001.  *   of twilight. eg, refraction to a level horizon is often taken to be about
  2002.  *   32/2(sun's semi-diameter) + 34(nominal refraction) = 50 arc minutes.
  2003.  *   better is:
  2004.  *      unrefract (pre, temp, 0.0, &a);    /* true alt of upper limb 
  2005.  *      a -= SUN_DIAM;            /* true alt of sun center
  2006.  * also used to find astro twilight by calling with dis of 18 degrees.
  2007.  * status:
  2008.  *   2: error
  2009.  *   1: never rises
  2010.  *   0: normal
  2011.  *  -1: circumpolar
  2012.  */
  2013. sunrs (mjd, lat, lng, dis, utcr, utcs, azr, azs, status)
  2014. float mjd, lat, lng, dis;
  2015. float *utcr, *utcs;
  2016. float *azr, *azs;
  2017. int *status;
  2018. {
  2019.     float lstr, lsts;    /* local sidereal times of rising/setting */
  2020.     float al = radhr(lng);    /* time correction for longitude */
  2021.     float tmp1, tmp2;    /* tmp */
  2022.     float mjd0, t;
  2023.  
  2024.     mjd0 = floor(mjd+0.5)-0.5; /* insure mjd0 is greenwich start of day */
  2025.  
  2026.     /* first approximation is to find rise/set times of a fixed object
  2027.      * in the position of the sun at local midday. the sun is not
  2028.      * fixed but moves about a degree per day so these are then refined.
  2029.      */
  2030.     fixedsunriset(mjd0+(12.-al)/24.,lat,dis,&lstr,&lsts,&tmp1,&tmp2,status);
  2031.     if (*status != 0) return;
  2032.  
  2033.     /* find a better approximation to the rising circumstances based on a
  2034.      * fixed object at the suns location at the first approximation of the
  2035.      * rise time.
  2036.      * N.B. more iterations are less than sp float precision.
  2037.      */
  2038.     gst_utc (mjd0, lstr, &t);
  2039.     fixedsunriset(mjd0+(t-al)/24.,lat,dis,&lstr,&tmp1,azr,&tmp2,status);
  2040.     if (*status != 0) return;
  2041.     gst_utc (mjd0, lstr, utcr);
  2042.     *utcr -= al*.9972696;    /* allow for sideral shift */
  2043.  
  2044.     /* find a better approximation to the setting circumstances based on a
  2045.      * fixed object at the suns location at the first approximation of the
  2046.      * setting time.
  2047.      */
  2048.     gst_utc (mjd0, lsts, &t);
  2049.     fixedsunriset (mjd0+(t-al)/24.,lat,dis,&tmp1,&lsts,&tmp2,azs,status);
  2050.     if (*status != 0) return;
  2051.     gst_utc (mjd0, lsts, utcs);
  2052.     *utcs -= al*.9972696;    /* allow for sideral shift */
  2053. }
  2054.  
  2055. /* find the local rise/set sidereal times and azimuths of an object fixed
  2056.  * at the ra/dec of the sun on the given mjd time as seen from lat.
  2057.  */
  2058. static
  2059. fixedsunriset (mjd, lat, dis, lstr, lsts, azr, azs, status)
  2060. float mjd, lat, dis;
  2061. float *lstr, *lsts;
  2062. float *azr, *azs;
  2063. int *status;
  2064. {
  2065.     float lsn, rsn;
  2066.     float deps, dpsi;
  2067.     float ra, dec;
  2068.  
  2069.     /* find ecliptic position of sun at mjd */
  2070.     sun (mjd, &lsn, &rsn);
  2071.  
  2072.     /* allow for 20.4" aberation and nutation */
  2073.     nutation (mjd, &deps, &dpsi);
  2074.         lsn += dpsi - degrad(20.4/3600.0);
  2075.  
  2076.     /* convert ecliptic to equatorial coords */
  2077.     ecl_eq (mjd, 0.0, lsn, &ra, &dec);
  2078.  
  2079.     /* find circumstances for given horizon displacement */
  2080.     riset (ra, dec, lat, dis, lstr, lsts, azr, azs, status);
  2081. }
  2082. xXx
  2083. echo extracting utc_gst.c
  2084. cat > utc_gst.c << 'xXx'
  2085. /* given a modified julian date, mjd, and a universally coordinated time, utc,
  2086.  * return greenwich mean siderial time, *gst.
  2087.  */
  2088. utc_gst (mjd, utc, gst)
  2089. float mjd;
  2090. float utc;
  2091. float *gst;
  2092. {
  2093.     float tnaught();
  2094.     static float lastmjd;
  2095.     static float t0;
  2096.  
  2097.     if (mjd != lastmjd) {
  2098.         t0 = tnaught (mjd);
  2099.         lastmjd = mjd;
  2100.     }
  2101.     *gst = 1.002737908*utc + t0;
  2102.     range (gst, 24.0);
  2103. }
  2104.  
  2105. /* given a modified julian date, mjd, and a greenwich mean siderial time, gst,
  2106.  * return universally coordinated time, *utc.
  2107.  */
  2108. gst_utc (mjd, gst, utc)
  2109. float mjd;
  2110. float gst;
  2111. float *utc;
  2112. {
  2113.     float tnaught();
  2114.     static float lastmjd;
  2115.     static float t0;
  2116.  
  2117.     if (mjd != lastmjd) {
  2118.         t0 = tnaught (mjd);
  2119.         range (&t0, 24.0);
  2120.         lastmjd = mjd;
  2121.     }
  2122.     *utc = gst - t0;
  2123.     range (utc, 24.0);
  2124.     *utc *= .9972695677;
  2125. }
  2126.  
  2127. static float
  2128. tnaught (mjd)
  2129. float mjd;    /* julian days since 1900 jan 0.5 */
  2130. {
  2131.     float dmjd;
  2132.     int m, y;
  2133.     float d;
  2134.     float t, t0;
  2135.  
  2136.     mjd_cal (mjd, &m, &d, &y);
  2137.     cal_mjd (1, 0., y, &dmjd);
  2138.     t = dmjd/36525;
  2139.     t0 = 6.57098e-2 * (mjd - dmjd) - 
  2140.          (24 - (6.6460656 + (5.1262e-2 + (t * 2.581e-5))*t) -
  2141.            (2400 * (t - (((float)y - 1900)/100))));
  2142.     return (t0);
  2143. }
  2144. xXx
  2145.  
  2146.